function I=CalcPCAAxis(ModelData)
%CALCPCAAXIS Calculate the average PCA direction for each ROI in the model data.
%
%   This function calculates the average principal component analysis (PCA)
%   direction for each region of interest (ROI) in the provided model data.
%   The PCA direction is determined by averaging the primary PCA direction and
%   its negated direction across all patients. The result is a structure with
%   the main PCA component as the third column, the secondary component as the
%   second column, and the least variation as the first column.
%
%   Syntax
%   I = CalcPCAAxis(ModelData)
%
%   Inputs
%   ModelData : A structure array containing model data for multiple patients.
%               Each element of the array must have a 'ROIName' field
%               containing the name of the ROI and a 'CoeffPCA' field
%               containing the PCA coefficients for that ROI.
%
%   Outputs
%   I : A structure array with fields 'RoiName' and 'Axis'.
%       'RoiName' contains the name of each ROI, and 'Axis' contains a
%       three-element vector representing the average PCA direction for that ROI.
%       The vector is ordered as [X, Y, Z], where X is the least variation,
%       Y is the secondary component, and Z is the main PCA component.
%
%   Example
%   Calculate the average PCA direction for each ROI in the model data
%   ModelData = ...; % Assume ModelData is properly structured
%   PCAAxis = CalcPCAAxis(ModelData);
%
%   See also
%   pca : MATLAB's built-in PCA function
%
%   Notes
%   This function assumes that the PCA coefficients are already calculated
%   for each ROI in the model data. The averaging process is not performed
%   during the initial model data collection.
%
%   Copyright (c) 2024, Odense University Hospital

roiList=unique({ModelData.ROIName});
I=struct();
for iROI=1:length(roiList)
  index = strcmp({ModelData.ROIName},roiList{iROI});

  % eigenvectorsPCA is a 3x3 matrix containing the 3D eigenvectors per
  % ModelData table entry
  eigenvectorsPCA = [ModelData(index).CoeffPCA];
  
  % taking only the primary axis (i.e. first column of the 3x3 matrix 
  % in eigenvectorsPCA and skipping the other two. 
  % unitvectorPrimaryPCA will contain the primary axis of all ModelData
  % entries of the ROI of the current iteration (defined by index
  % variable). This means the dimension is 3xN where N is the number of
  % table entries.
  unitvectorPrimaryPCA = eigenvectorsPCA(:,1:3:end); 

  % unitvectorPrimaryNegatedPCA is relevant for allowing invariance along
  % an axis such that the + direction and - direction are not treated
  % differently in the averagedDirection computation
  unitvectorPrimaryNegatedPCA = [ ...
      unitvectorPrimaryPCA, ...
      -unitvectorPrimaryPCA ...
      ];
  
  % DEPRECATED:
  % The next line is basically equivalent to
  % averagedDirection = pca(unitvectorPrimaryNegatedPCA');
  % however, note that pca() does curiously enough not standardize the data
  % so that needs to be done actively (i.e. pca(normalize(myData))).
  % We use the manual approach, that is, calculating the eigenvectors of
  % the correlation matrix of the standardized data as this always returns
  % three vectors for 3D data. pca() gives less eigenvectors if the
  % system is underdetermined/undetermined.
  % [averagedDirection, ~] = eig(corrcoef(normalize(unitvectorPrimaryNegatedPCA')));
  % averageZ= averagedDirection(:,3);
  % averageY= averagedDirection(:,2);
  % averageX = cross(averageY,averageZ);
  
  averagedDirection = pca(unitvectorPrimaryNegatedPCA');
  averageZ = averagedDirection(:,1);
  averageY = averagedDirection(:,2);
  averageX = cross(averageY, averageZ);

  I(iROI).RoiName = roiList{iROI};
  I(iROI).Axis = [averageX, averageY, averageZ];
end
end