classdef Scapula < handle
    %SCAPULA Summary of this class goes here
    %   This class defines the scapula. Landmarks are used to define its
    %   coordinate system. It includes the glenoid object.

    properties
        angulusInferior % Landmark Angulus Inferior read from amira
        angulusInferiorLocal % Same as above, expressed the scapular coordinate system rather than the CT coordinate system
        trigonumSpinae % Landmark Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine read from amira
        processusCoracoideus % Landmark Processus Coracoideus (most lateral point) read from amira
        acromioClavicular % Landmark Acromio-clavicular joint (most lateral part on acromion)read from amira
        angulusAcromialis % Landmark Angulus Acromialis (most laterodorsal point) read from amira
        spinoGlenoidNotch % Landmark Spino-glenoid notch read from amira
        pillar % 5 landmarks on the pillar
        groove % 5 landmarks on the scapula (supraspinatus) groove

        coordSys % Scapular coordinate system
        plane % Plane class
        segmentation % either none 'N', manual 'M' or automatic 'A'

        surface % surface points and triangles of the scapula
        glenoid % Glenoid class
        acromion % Acromion class
        comment % any comment
    end

    properties (Hidden = true)
        SCase
    end

    methods (Access = ?Shoulder)  % Only Shoulder is allowed to construct a Scapula
        function obj = Scapula(SCase)
            %SCAPULA Construct an instance of this class
            %   Constructor of the scapula object, sets all properties to
            %   zero.
            obj.SCase = SCase;

            obj.angulusInferior       = [];
            obj.angulusInferiorLocal  = [];
            obj.trigonumSpinae        = [];
            obj.processusCoracoideus  = [];
            obj.acromioClavicular     = [];
            obj.angulusAcromialis     = [];
            obj.spinoGlenoidNotch     = [];
            obj.pillar                = [];
            obj.groove                = [];

            obj.coordSys              = CoordinateSystem;
            obj.plane                 = Plane();
            obj.segmentation          = 'N';

            obj.surface               = [];
            obj.glenoid               = Glenoid(SCase);
            obj.glenoid.scapula       = obj;
            obj.acromion              = Acromion(SCase);
            obj.acromion.scapula      = obj;

            obj.comment = '';
        end
    end

    methods
        function outputArg = load(obj)
            % LOAD Load segmented surface and landmnarks
            %   Load the segmented scapula surface (if exist) and the
            %   scapula landmarks (6 scapula, 5 groove, 5 pillar) from
            %   amira directory.

            SCase = obj.SCase;
            SCaseId4C = SCase.id4C;
            amiraDir = SCase.dataAmiraPath;
            matlabDir = SCase.dataMatlabPath;
            outputArg = 1; % Should report on loading result

            % Scapula (6) landmarks
            fileName = ['ScapulaLandmarks' SCaseId4C '.landmarkAscii'];
            fileName = [amiraDir '/' fileName];

            if exist(fileName,'file') == 2
                importedData = importdata(fileName, ' ', 14);

            else
                disp(['MATLAB:rmpath:DirNotFound',...
                    'Could not find file: ',fileName]);
                outputArg = 0;
                return;
            end
            % Check that imported data are well formed
            try
                landmarks = importedData.data;
            catch
                outputArg = 0;
                return;
            end

            obj.angulusInferior      = landmarks(1,:);
            obj.trigonumSpinae       = landmarks(2,:);
            obj.processusCoracoideus = landmarks(3,:);
            obj.acromioClavicular    = landmarks(4,:);
            obj.angulusAcromialis    = landmarks(5,:);
            obj.spinoGlenoidNotch    = landmarks(6,:);

            % Scapula suprapinatus groove (5) landmarks
            fileName = ['ScapulaGrooveLandmarks' SCaseId4C '.landmarkAscii'];
            fileName = [amiraDir '/' fileName];
            if exist(fileName,'file') == 2
                importedData = importdata(fileName, ' ', 14);
                % Check that imported data are well formed
                try
                    landmarks = importedData.data;
                catch
                    outputArg = 0;
                    return;
                end
            else
                disp(['MATLAB:rmpath:DirNotFound',...
                    'Could not find file: ',fileName]);
                outputArg = 0;
            end
            % Check groove landmarks validity
            if length(landmarks) > 5
                landmarks = landmarks(1:5,:);
                warning(['More than 5 groove landmarks(' SCase.id ')']);
            elseif length(landmarks) < 5
                error('Less than 5 groove landmarks');
            end

            obj.groove = landmarks;

            % Scapula pillar (5) landmarks
            fileName = ['ScapulaPillarLandmarks' SCaseId4C '.landmarkAscii'];
            fileName = [amiraDir '/' fileName];
            if exist(fileName,'file') == 2
                importedData = importdata(fileName, ' ', 14);
                % Check that imported data are well formed
                try
                    landmarks = importedData.data;
                catch
                    outputArg = 0;
                    return;
                end
            else
                disp(['MATLAB:rmpath:DirNotFound',...
                    'Could not find file: ',fileName]);
                outputArg = 0;
            end
            obj.pillar = landmarks;

            % Import scapula surface
            % Try manual segmentation in amira folder
            fileName = ['scapula_' SCaseId4C '.stl'];
            fileName = [amiraDir '/' fileName];
            if exist(fileName,'file') == 2
                try
                    [points,faces,~]=loadStl(fileName,1);
                    obj.surface.points = points;
                    obj.surface.faces = faces;
                    obj.segmentation = 'M'; % Manual segmentation from amira
                catch
                    obj.segmentation = 'E'; % Error in loading
                end
            else
                obj.segmentation = 'N'; % No segmentation
            end

            obj.measurePlane;
        end

        function outputArg = loadAuto(obj)
            % Load 2 files obtained by SSM auto segmentation
            % scapulaSurfaceAuto{L/R}.ply
            % scapulaLandmarksAuto{L/R}.mat

            outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise

            SCase = obj.SCase;
            SCaseId4C = SCase.id4C;
            matlabDir = SCase.dataMatlabPath;

            side = SCase.shoulder.side; % We might no have it yet !!

            if isempty(side)
              side = 'R';
            end

            % Try loading auto segmentation from matlab dir
            fileName = ['scapulaSurfaceAuto' side '.ply'];
            fileName = [matlabDir '/' fileName];

            if ~exist(fileName,'file')
              side = 'L';
              fileName = ['scapulaSurfaceAuto' side '.ply'];
              fileName = [matlabDir '/' fileName];
            end


            %side = SCase.shoulder.side;

            if exist(fileName,'file') == 2
                try
                    face_index_start = 1;
                    [points,faces] = loadPly(fileName,face_index_start); % % Load ply file

                    % ptCloud = pcread(fileName); % load the pointCloud object in the ply file
                    % points = ptCloud.Location; % get the array of (x,y,z) points

                    obj.surface.points = points;
                    obj.surface.faces = faces;
                    obj.segmentation = 'A'; % Auto segmentation from matlab
                catch
                    obj.segmentation = 'E'; % Error on loading
                    outputArg = 0;
                end
            else
                obj.segmentation = 'N'; % No segmentation file
                outputArg = 0;
            end

            % Try loading auto scapula landmarks from matlab dir
            filename = ['scapulaLandmarksAuto' side '.mat'];
            filename = [matlabDir '/' filename];
            if exist(filename,'file') == 2
                try
                    % Load mat file with scapula landmarks
                    load(filename, 'ScapulaLandmarks');
                    % Set loaded landmarks to scapulaAuto (obj)
                    obj.angulusInferior       = ScapulaLandmarks.angulusInferior;
                    obj.trigonumSpinae        = ScapulaLandmarks.trigonumSpinae;
                    obj.processusCoracoideus  = ScapulaLandmarks.processusCoracoideus;
                    obj.acromioClavicular     = ScapulaLandmarks.acromioClavicular;
                    obj.angulusAcromialis     = ScapulaLandmarks.angulusAcromialis;
                    obj.spinoGlenoidNotch     = ScapulaLandmarks.spinoGlenoidNotch;
                    obj.pillar                = ScapulaLandmarks.pillar;
                    obj.groove                = ScapulaLandmarks.groove;
                catch
                    disp(SCaseId4C); % For debug (2 cases)
                    outputArg = 0;
                end
            else
                outputArg = 0;
            end

            obj.measurePlane;
        end

        function output = coordSysSet(obj)
          % Choose your method with the next command
          obj.coordSys.setSystemWithScapulaLandmarks(obj);
          % obj.coordSys.setSystemWith;

          % To enable retro compatibility with scapula.coordSysSet former
          % implementation, the two following commands are called here
          obj.setShoulderSideWithCoordinateSystem;
          obj.angulusInferiorLocal = obj.coordSys.express(obj.angulusInferior);
          output = 1;
        end

        function measurePlane(obj)
          % Scapular plane is fitted on 3 points (angulusInferior,
          % trigonumSpinae, most laretal scapular goove landmark).
          inferior = obj.angulusInferior;
          medial = obj.trigonumSpinae;
          mostLateralGrooveIndex = findLongest3DVector(medial-obj.groove);
          mostLateralGroovePoint = obj.groove(mostLateralGrooveIndex(1),:);
          obj.plane.fit([inferior; medial; mostLateralGroovePoint]);

          anterior = obj.processusCoracoideus;
          posterior = obj.angulusAcromialis;
          obj.plane.normal = orientVectorToward(obj.plane.normal,(anterior-posterior));
        end

        function setShoulderSideWithCoordinateSystem(obj)
          % The scapula coordinate system is right-handed for the right shoulder only
          if obj.coordSys.isRightHanded
            obj.SCase.shoulder.side = 'R';
            return
          elseif obj.coordSys.isLeftHanded
            obj.SCase.shoulder.side = 'L';
            return
          end
          warning(['Couldn''t find the shoulder side by evaluating the scapula',...
                   'coordinate system.']);
        end

    end
end