classdef Acromion < handle
    %ACROMION Properties and methods associted to the acromion.
    %   Detailed explanation goes here
    
    properties
        AI   % Acromion Inddex (doi:10.2106/JBJS.D.03042)
        CSA  % Critital Shoulder Angle (doi:10.1302/0301-620X.95B7.31028)
        PSA  % Posterior Slope Angle (10.1016/S1058-2746(05)80036-9)
        PSL  % Length of segment beteeen AA and AC
        AAA  % Angulus Angle Angle (Angle between AA-origin and PA axis)
        AAL  % Length of segment between origine and AA
        AAx  % PA (x) position of AA
        AAy  % IS (y) position of AA
        AAz  % ML (z) position of AA
        ACx  % PA (x) position of AC
        ACy  % IS (y) position of AC
        ACz  % ML (z) position of AC
        
        comment
    end
    
    properties (Hidden = true)
        SCase
    end
    
    methods (Access = ?Scapula)  % Only Scapula is allowed to construct a Acromion
        function obj = Acromion(SCase)
            %ACROMION Construct an instance of this class
            %   Instance of the acromion objet set all properties to zero;
            obj.AI   = [];
            obj.CSA  = [];
            obj.PSA  = [];
            obj.PSL  = [];
            obj.AAA  = [];
            obj.AAL  = [];
            obj.AAx  = [];
            obj.AAy  = [];
            obj.AAz  = [];
            obj.ACx  = [];
            obj.ACy  = [];
            obj.ACz  = [];
            obj.comment = '';

            obj.SCase = SCase;
        end
    end
    
    methods
        function outputArg = morphology(obj)
            %MORPHOLOGY Performs morphology analysis of the acromion.
            % It caluculates acromion index (AI), critical shoulder angle
            % (CSA), and posterior slope (PS). The glenoid and scapula
            % surface are re-oriented (rotated and translated) in the
            % scapular coordinate system. For AI nd CSA, glenoid and
            % scapula points are projected in the scapular plane, which is
            % [1 0 0] after the re-orientation.
            % TODO:
            % Might be removed anatomy for concistency

            SCase = obj.SCase;
            
            %% Get local coordinate system from parent scapula
            origin = SCase.shoulder.scapula.coordSys.origin;
            xAxis  = SCase.shoulder.scapula.coordSys.PA;
            yAxis  = SCase.shoulder.scapula.coordSys.IS;
            zAxis  = SCase.shoulder.scapula.coordSys.ML;
            
            % Rotation matrix to align to scapula coordinate system
            R = [xAxis' yAxis', zAxis'];
            
            %% Scapula, glenoid, and scapula landmarks in scapula coordinate system
            
            % Scapula surface alignment is done in next section, for the
            % caulation of AI, because of mix between auto, manual and no
            % segmentation. This should be corrected when there will be two
            % scapula object per shoulder (manual and auto).
            
            % Glenoid surface in scapular coordinate system
            glenSurf = (SCase.shoulder.scapula.glenoid.surface.points - origin) * R;
            % Glenoid center in scapular coordinate system
            glenCenter = (SCase.shoulder.scapula.glenoid.center - origin) * R;
            % Scapula landmarks in scapular coordinate system
            AC = (SCase.shoulder.scapula.acromioClavicular - origin) * R; % Acromio-clavicular landmark
            AA = (SCase.shoulder.scapula.angulusAcromialis - origin) * R; % Angulus acromialis landmark
            obj.ACx = AC(1);
            obj.ACy = AC(2);
            obj.ACz = AC(3);
            obj.AAx = AA(1);
            obj.AAy = AA(2);
            obj.AAz = AA(3);
                         
            %% Acromion Index (AI)
            % Adapted by AT (2018-07-18) from Valerie Mafroy Camine (2017)
            
            % AI = GA/GH, where:
            % GA: GlenoAcromial distance in scapula plane = distance from
            % AL to glenoid principal axis
            % GH: GlenoHumeral distance = 2*HHRadius, humeral head diameter (radius*2)
            % AL: most lateral point of the acromion
            
            % Get all required data, aligned in scapular plane, and
            % project in scapular plane.
            
            ScapPlaneNormal = [1 0 0]; % Normal of the scapular plane in the scapular coordinate system
            PlaneMean = [0 0 0];       % Origin of scapular system in scapular coordinate system
            
            % Project glenoid surface in scapular plane
            glenSurf = project2Plane(glenSurf,ScapPlaneNormal,PlaneMean,size(glenSurf,1));
            
            % Project glenoid center in scapular plane
            glenCenter = project2Plane(glenCenter,ScapPlaneNormal,PlaneMean,size(glenCenter,1));
            
            % If scapula is segmented, get AL from most lateral part of the
            % scapula, otherwise use acromio-clavicular landmark
            % Get segmentation propertiy from parent scapula
            segmentedScapula = SCase.shoulder.scapula.segmentation;
            segmentedScapula = strcmp(segmentedScapula,'A') || strcmp(segmentedScapula,'M');
            if segmentedScapula == 1
                % Rotate scapula surface to align the scapular plane with
                % the YZ plane, and then take the max Z (lateral) point
                % Transform scapula surface to scapula coordinate system
                scapSurf = SCase.shoulder.scapula.surface.points;
                scapSurf = (scapSurf - origin) * R;
                scapSurf = project2Plane(scapSurf,ScapPlaneNormal,PlaneMean,size(scapSurf,1));
                
                % Take the most lateral point of the scapula, assuming that
                % this point is the most lateral point of the acromion
                [~,ALpositionInScapula] = max(scapSurf(:,3));
                AL = scapSurf(ALpositionInScapula,:);
            else
                % No scapula points, we approximate AL with the acromioClavicular
                % landmark, in the scapula coordinate system
                AL = AC; % Acroomio-clavicalar scapula landmark
                AL = project2Plane(AL,ScapPlaneNormal,PlaneMean,size(AL,1));
            end
            
            % Find glenoid principal axis with most superior and most inferior projected
            % glenoid points
            % AT: This method is not ideal. PCA would be better. It is also
            % used below by the CSA
            glenPrinAxis = [...
                glenSurf(glenSurf(:,2) == min(glenSurf(:,2)),:) ;...
                glenSurf(glenSurf(:,2) == max(glenSurf(:,2)),:)...
                ];
            
            % Compute GA (distance from AL to glenoid principal axis)
            % Most inferior point of the scapula surface
            IG = glenPrinAxis(1, :);
            % Most superior point of the scapula surface
            SG = glenPrinAxis(2, :);
            GA = norm(cross(SG - IG, AL - IG))/norm(SG - IG);
            
            % GH (Humeral head diameter)
            % get humeral head radius from associated humerus
            GH = 2 * SCase.shoulder.humerus.radius;
            
            % Acromion index
            obj.AI = GA/GH;
            
            %% Critical Shoulder Angle (CSA)
            % Adapted by AT (2018-07-18) from Bharath Narayanan (2018)
            % [radCSA, degCSA] = computeCSA(GlenoidPtsinScapulaPlane, ALinScapulaPlane);
            
            % Vectors connecting IG to SG, and IG to AL
            IGSG = SG - IG;
            IGAL = AL - IG;
            CSA = vrrotvec(IGSG,IGAL);
            CSA = CSA(4);
            CSA = rad2deg(CSA);
            
            obj.CSA = CSA;
            
            %% Posterior Slope Angle and length (PSA & PSL)
            % By AT (2018-08-10)
            
            % Project AA & AC in PA-IS plane
            AAxy = [AA(1:2) 0]; % Juste take x and y component
            ACxy = [AC(1:2) 0]; % Juste take x and y component
            PSv  = ACxy - AAxy; % Posterior slope vector
            ISv  = [1 0 0]; % IS axis
            PSA = vrrotvec(PSv, ISv);
            PSA = PSA(4);
            PSA = rad2deg(PSA);
            
            PSL = norm(PSv);
            
            obj.PSA = PSA;
            obj.PSL = PSL;
            
            %%  Acromial Angle Angle and length (AAA & AAL)
            % By AT (2018-09-13)
            
            % Vector between scapular origin and AA in the plane PA-IS
            AAv = [AA(1:2) 0];
            
            % Angle between AAvect and PA axis
            PAv = [1 0 0]; % PS axis
            AAA = vrrotvec(AAv, PAv);
            AAA = AAA(4);
            AAA = rad2deg(AAA);
            AAA = 180 - AAA;
            AAL = norm(AAv);
            
            obj.AAA = AAA;
            obj.AAL = AAL;
            
            %% We might save an image of AI, CSA, PS
            
            outputArg = 1;
        end
    end
end