classdef (Abstract) Shoulder < handle

  properties
    side % R or L
    scapula
    humerus
    muscles
    CTscan
    comment
  end



  properties (Hidden = true)
    SCase
  end



  methods (Abstract, Access = protected)
    createScapula(obj);
  end



  methods
    function obj = Shoulder(SCase)
      obj.side    = '';
      obj.humerus = Humerus(obj);
      obj.muscles = Muscles(obj);
      obj.CTscan  = '';
      obj.comment = '';

      obj.SCase   = SCase;

      obj.createScapula;
    end

    function output = plot(obj)
      % Create a figure and plot inside:
      % the scapula surface, coordSys, groove points and landmarks;
      % the glenoid surface, sphere and centerLine;
      % the humerus sphere and centerLine.

      figure;
      hold on;


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %% Scapula
      % Plot the segmented surface
      if not(isempty(obj.scapula.surface))
          points = obj.scapula.surface.points;
          faces  = obj.scapula.surface.faces;
          x = points(:,1);
          y = points(:,2);
          z = points(:,3);
          trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud');
      end

      % Plot scapula groove
      if not(isempty(obj.scapula.groove))
        points = obj.scapula.groove;
        x = points(:,1);
        y = points(:,2);
        z = points(:,3);
        scatter3(x,y,z,100,'blue');
      end

      % Plot scapula markers
      points = [...
          obj.scapula.angulusInferior;
          obj.scapula.trigonumSpinae;
          obj.scapula.processusCoracoideus;
          obj.scapula.acromioClavicular;
          obj.scapula.angulusAcromialis;
          obj.scapula.spinoGlenoidNotch...
          ];
      x = points(:,1);
      y = points(:,2);
      z = points(:,3);
      scatter3(x,y,z,100,'blue','LineWidth',2);


      % Plot scapular plane
      plot(obj.scapula.plane,100);

      % Plot scapular axis
      pt1 = obj.scapula.coordSys.origin;
      axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]);
      pt2 = pt1 - obj.scapula.coordSys.ML * axisLength;
      x = [pt1(1), pt2(1)];
      y = [pt1(2), pt2(2)];
      z = [pt1(3), pt2(3)];
      plot3(x,y,z,'Color','blue','LineWidth', 5);

      % Plot scapular coordinate system
      axisLength = 50; % Length of the coord. sys. axes
      pt0 = obj.scapula.coordSys.origin;
      pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X
      pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y
      pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z
      x = [pt0(1), pt1(1)];
      y = [pt0(2), pt1(2)];
      z = [pt0(3), pt1(3)];
      plot3(x,y,z,'Color','red','LineWidth', 5);
      x = [pt0(1), pt2(1)];
      y = [pt0(2), pt2(2)];
      z = [pt0(3), pt2(3)];
      plot3(x,y,z,'Color','green','LineWidth', 5);
      x = [pt0(1), pt3(1)];
      y = [pt0(2), pt3(2)];
      z = [pt0(3), pt3(3)];
      plot3(x,y,z,'Color','blue','LineWidth', 5);

      % plot scapular axis on glenoid center
      axisLength = 50;
      pt1 = obj.scapula.glenoid.center;
      pt2 = pt1 + obj.scapula.coordSys.ML*axisLength;
      x = [pt1(1), pt2(1)];
      y = [pt1(2), pt2(2)];
      z = [pt1(3), pt2(3)];
      plot3(x,y,z,'Color','magenta','LineWidth', 5);


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Glenoid

      % Plot the glenoid surface
      % We might move them in the direction of th glenoid center
      % instead of the ML axis

      points = obj.scapula.glenoid.surface.points +...
          obj.scapula.coordSys.ML * 0.2;
      faces  = obj.scapula.glenoid.surface.faces;
      x = points(:,1);
      y = points(:,2);
      z = points(:,3);
      trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none');


      % plot glenoid centerline
      pt1 = obj.scapula.glenoid.center;
      pt2 = pt1 + obj.scapula.glenoid.centerLine;
      x = [pt1(1), pt2(1)];
      y = [pt1(2), pt2(2)];
      z = [pt1(3), pt2(3)];
      plot3(x,y,z,'Color','magenta','LineWidth', 5);

      % We might add spherical cap

      % plot glenoid sphere
      n = 20; % number of lines
      [X,Y,Z] = sphere(n);
      X = X*obj.scapula.glenoid.radius;
      Y = Y*obj.scapula.glenoid.radius;
      Z = Z*obj.scapula.glenoid.radius;
      Xc = obj.scapula.glenoid.center(1);
      Yc = obj.scapula.glenoid.center(2);
      Zc = obj.scapula.glenoid.center(3);
      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Humerus

      % plot humeral head sphere
      n = 20; % number of lines
      [X,Y,Z] = sphere(n);
      X = X*obj.humerus.radius;
      Y = Y*obj.humerus.radius;
      Z = Z*obj.humerus.radius;
      Xc = obj.humerus.center(1);
      Yc = obj.humerus.center(2);
      Zc = obj.humerus.center(3);
      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');


      % plot humeral head centerline
      pt1 = obj.scapula.glenoid.center;
      pt2 = obj.humerus.center;
      x = [pt1(1), pt2(1)];
      y = [pt1(2), pt2(2)];
      z = [pt1(3), pt2(3)];
      plot3(x,y,z,'Color','magenta','LineWidth', 5);


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Graph properties

      axis off; % remove axis
      lightPosition = obj.scapula.glenoid.center + ...
          (...
          obj.scapula.coordSys.ML + ...
          obj.scapula.coordSys.IS + ...
          obj.scapula.coordSys.PA ...
          ) * 100;
      light('Position',lightPosition,'Style','local');
      grid on;
      xlabel('X (CT)');
      ylabel('Y (CT)');
      zlabel('Z (CT)');
      axis square;
      axis vis3d; % fixed limits and scaling
      viewPoint = obj.scapula.coordSys.ML + ...
                  obj.scapula.coordSys.IS + ...
                  obj.scapula.coordSys.PA;
      viewPoint = obj.scapula.coordSys.PA;
      view(viewPoint); % view from lateral side
      % Align IS axis with window vertical axis
      camup(obj.scapula.coordSys.IS);

      figureHandle = gca;
      figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis
      figureHandle.Box = 'On';
      set(gcf,'color','w'); % Set background color to white

      output = 1;
    end
  end
end