function I=CalcBasicFeatures(CT,Mask,varargin)
% This function calculates some of the basic features for a given Mask. 
% Image data (named CT) is used to define the coordinate system of the Mask (voxel coordinates).
% Mask is output from StructMask_OUH
% The features that are calculated are:
% 1)	Volume
% 2)	Surface area
% 3)	Compactness
% 4)	Center of mass
% 5)	Longest ROI distance
% 6)	Length of PCA vectors
% 7)	Direction angles of the largest PCA component (orientation in space)
% 8)	Distance from centre of mass to surface (negative if the centre of mass is outside Mask)
% All values are returned in the struct I
%
% Example Usage:
% files=SortDicomFiles_OUH(dataPath);
% MR=LoadDicomData_OUH(files,5,1);
% RTStruct=LoadDicomData_OUH(files,4,1);
% mask=StructMask_OUH(RTStruct,MR,iRois);
% basicFeatures = CalcBasicFeatures(MR,mask);

% Input parser setup
inputParam = inputParser;
addRequired(inputParam, 'CT', @(x) iscell(x));
addRequired(inputParam, 'Mask', @(x) iscell(x));

addParameter(inputParam, 'fastMode', false, @(x) islogical(x) || isnumeric(x))

% Parse inputs
parse(inputParam, CT, Mask, varargin{:});
inputValues = inputParam.Results;


[x,y,z] = GetVoxelPositions_OUH(CT);
% Number of Voxels in Mask
NrVoxels = sum(Mask{2}(:));

% Calculate Volume (mm3)
% Does so by calculating the volume of a single voxel
% then multiplying it with the number of voxels in the ROI Mask
voxelSize = (x(2)-x(1))*(y(2)-y(1))*(z(2)-z(1));
I.Volume = NrVoxels*voxelSize;

    if I.Volume > 0
        % Calculate the center of mass
        [mesh1,mesh2,mesh3] = meshgrid(x,y,z);
        index = find(Mask{2}==1); % get all indices where binary mask of ROI is 1 
        positions = [mesh1(index), mesh2(index), mesh3(index)]; % filter mesh to only the meshes where binary mask is 1
        meanPositions = mean(positions);
        I.xcm = meanPositions(1); % x-pos of center of mass
        I.ycm = meanPositions(2); % y-pos of center of mass
        I.zcm = meanPositions(3); % z-pos of center of mass
        
        % PCA Explanation
        % Principal component analysis can be performed in the following way:
        % 1)	The initial matrix is named positions
        % 2)	Subtract the mean from each column, name that AMatrix
        % 3)	Compute the correlation matrix, name that CMatrix
        % 4)	Calculate the eigenvectors and eigenvalues
        % 5)	Sort to have eigenvectors ordered after descending eigenvalues
        % 6)	Change direction of eigenvector such that the third component is always positive
        % The above can be inplemented as shown in comments below or by using
        % matlab pca as being done below.
        % temp=mean(positions);
        % Amatrix=positions-repmat(temp,size(positions,1),1);
        % Cmatrix=(Amatrix'*Amatrix)/(size(positions,1)-1);
        % [v,d]=eig(Cmatrix);
        % temp=diag(d);
        % [~,ix]=sort(temp,'descend');
        % v=v(:,ix);
        % d=temp(ix);
        % signval=sign(v(end,:));
        % v=v.*repmat(signval,3,1);
        % After the next three lines coeff is equal to v (except for sign) and latent=d
        
        % Using Matlab PCA algorithm to rotate
        % the coordinate system of the datapoints
        % that make up a ROI defined in an RTStruct
        % (i.e. DICOM tag (3006, 0050) Data Points).
        % This will give us three eigenvectors from the PCA 
        % (i.e. that axis of the ROI) and three eigenvalues
        % from the PCA (i.e. the length of the ROI along the axis).
        % We can use the axis furthermore to calculate the euler angles.
        % This allows us to compare angles between the same ROI in different
        % patients (not part of this code).
        % Score returns the datapoints (i.e. position) after rotation on
        % the new PCA axis.
        [coeff,score,latent] = pca(positions);
        I.lengthPCA1 = latent(1);
        I.lengthPCA2 = latent(2);
        I.lengthPCA3 = latent(3);
        I.transformedCoord = score;
        
        % If an ROI is only 2D PCA3 might be zero, in that case, replace the value
        % with the smallest dimension of the voxels
        if I.lengthPCA3<(min([x(2)-x(1),y(2)-y(1),z(2)-z(1)])/10)
          I.lengthPCA3=min([x(2)-x(1),y(2)-y(1),z(2)-z(1)]);
        end  
        I.CoeffPCA=coeff; % NOTE: By default Matlab's PCA does not output coeff values that span a right handed coordinate system. (because right-handedness is not defined for N-dimensional problems only for 3-Dimensional)
        
        % Calculate the surface of ROI
        % The area of a triangle can be calculated as half the magnitude of the cross-product of two of its sides
        % We use isosruface for this which uses the marching cubes
        % algorithm and outputs a trangulated mesh that approximates the
        % isosurface at level 0.5 of a binary volume. The pros of this is
        % that acccuracy is super high  while the con is the computational
        % high cost.
        % ORIGINAL:
        if ~inputValues.fastMode
          FV = isosurface(x,y,z, Mask{2}, 0.5);
          V = FV.vertices;
          F = FV.faces;
        else
          % Faster method by downsampling the mask
          % Specify downsampling factor (e.g. 2 or 4)
          ds = 2;
  
          % Downsample mask for speed (reduce resolution)
          maskDs = Mask{2}(1:ds:end, 1:ds:end, 1:ds:end);
  
          % Create spatial coordinates
          xDs = x(1:ds:end);
          yDs = y(1:ds:end);
          zDs = z(1:ds:end);
  
          % Use isosurface
          FV = isosurface(xDs, yDs, zDs, maskDs, 0.5);
          V = FV.vertices;
          F = FV.faces;
        end

        % faster alternative but not as accurate with positioning
        % [F,V] = extractIsosurface(Mask{2}, 0.5);

        A =(V(F(:, 2), :) - V(F(:, 1), :));
        B =(V(F(:, 3), :) - V(F(:, 1), :));
        C = cross(A, B, 2);
        I.Area = 0.5 * sum(sqrt(sum(C.^2, 2)));
                
        % Calculate the number of islands in a structure
        % SUGGESTION: This could be potentially used for 
        % area calculation and CM to surface (not sure about the latter).
        % It would speed up the calculation substantially as regionprops3()
        % is much faster than isosurfaces()
        % Uses as default connectivity of 26 (see
        % https://se.mathworks.com/help/images/pixel-connectivity.html).
        % If we want to use another connectivity (i.e. only face touches
        % are connected which is equivalent to connectivity = 6) then we
        % need to use bwconncomp()
        geometricProperties = regionprops3(Mask{2}>0, "all");
        

        % Solidity is the ratio between a convex hull around the structure
        % and the actual structure within. Could be potentially added.
        I.Solidity = mean(geometricProperties.Solidity);

        % each connected component is a row in the regionprops3 table. Thus for each disconnected structure an additional row is added
        % The definition should be that if there is just one main
        % structure, the island parameter is 0. If there are two structures
        % (i.e. one main component and one island) the island value should
        % be 1.
        I.Islands = height(geometricProperties) - 1;
        
        % Compactness measures how closely packed an object is. 
        % A more compact object has a smaller radius of gyration.
        % For example:
        % A sphere is the most compact object.
        % A cube is less compact than a sphere because of its corners and edges.
        I.Compactness=36*pi*I.Volume^2/I.Area^3; 
        
        % Distance form center of mass to surface
        temp=V-repmat([I.xcm,I.ycm,I.zcm],size(V,1),1);
        distsurface=sqrt(min(sum(temp.*temp,2)));
        indexx=find(abs(x-I.xcm)==min(abs(x-I.xcm))); indexx=indexx(1);
        indexy=find(abs(y-I.ycm)==min(abs(y-I.ycm))); indexy=indexy(1);
        indexz=find(abs(z-I.zcm)==min(abs(z-I.zcm))); indexz=indexz(1);
        if Mask{2}(indexy,indexx,indexz)==0
          distsurface=-distsurface;
        end
        I.CMToSurface=distsurface;
        
        % % This used to calculate the maximum distance between any two
        % % arbitrary points.
        % distances = pdist2(V, V, 'euclidean');
        % maxDistance = max(distances(:));
        % I.maxLength = maxDistance;
        
        % Maximum distance from center of mass to a datapoint
        % This used to calculate the maximum distance between any two
        % arbitrary points. For efficency we changed it to calculate the
        % maximum distance between any vertice with the center of mass.
        % This should allow to find spots far away from the ROI.
        distVerticesCenterOfMass = V - [I.xcm, I.ycm, I.zcm];
        maxLength=max(sum(distVerticesCenterOfMass.*distVerticesCenterOfMass,2));
        I.MaxLength=sqrt(maxLength);
        
        % intensityMask = CT{1,2} .* uint16(Mask{1,2});
        % minimumIntensity = min(intensityMask(intensityMask > 0), [], "all")
        % maximumIntensity = max(intensityMask(intensityMask > 0), [], "all")
        % meanIntensity = mean(intensityMask(intensityMask > 0), "all")
        % medianIntensity = median(intensityMask(intensityMask > 0), "all")

    else
        I = NaN;
        disp('ROI has Volume 0. Skipping further feature calculations...')
    end
end