classdef RotatorCuff < handle

  properties
    SC
    SS
    IS
    TM
    anteroPosteriorImbalance = nan;

    shoulder;
  end

  methods
    function obj = RotatorCuff(shoulder)
      obj.shoulder = shoulder;

      obj.SC = Muscle(obj, "SC");
      obj.SS = Muscle(obj, "SS");
      obj.IS = Muscle(obj, "IS");
      obj.TM = Muscle(obj, "TM");

      [~, ~] = mkdir(obj.dataPath);
      [~, ~] = mkdir(obj.dataSlicesPath);
    end

    function output = dataPath(obj)
      output = fullfile(obj.shoulder.dataPath, "muscles");
    end

    function output = dataSlicesPath(obj)
      output = fullfile(obj.dataPath, "oblique_slices");
    end

    function output = summary(obj)
      summary = [];
      summary = [summary; obj.SC.summary];
      summary = [summary; obj.SS.summary];
      summary = [summary; obj.IS.summary];
      summary = [summary; obj.TM.summary];
      output = summary;
    end

    function sliceAndSegment(obj, varargin)
      obj.createSliceMatthieu();
      obj.segmentMuscles("rotatorCuffMatthieu", "autoMatthieu");
    end

    function output = measure(obj, sliceName, segmentationSet)
      try
        obj.measureDegenerations(sliceName, segmentationSet);
        obj.measureInsertions();
        obj.measureAnteroPosteriorImbalance();
      end
      output = obj.summary;
    end

    function measureDegenerations(obj, sliceName, segmentationSet)
      for muscleName = ["SC", "SS", "IS", "TM"]
        obj.(muscleName).measureDegeneration(sliceName, segmentationSet);
      end
    end

    function measureInsertions(obj)
      humerus = obj.shoulder.humerus;
      humeralHead = Sphere(humerus.center, humerus.radius);

      % IS humeral insertion estimation
      obj.IS.measureInsertionTangentToSphere(...
        humeralHead,...
        obj.shoulder.scapula.coordSys.IS,...
        -obj.shoulder.scapula.coordSys.PA);
      obj.IS.measureForceVector();

      % TM humeral insertion estimation
      obj.TM.measureInsertionTangentToSphere(...
        humeralHead,...
        obj.shoulder.scapula.coordSys.IS,...
        -obj.shoulder.scapula.coordSys.PA);
      obj.TM.measureForceVector();

      % SC humeral insertion estimation
      obj.SC.measureInsertionTangentToSphere(...
        humeralHead,...
        obj.shoulder.scapula.coordSys.IS,...
        obj.shoulder.scapula.coordSys.PA);
      obj.SC.measureForceVector();
    end

    function measureAnteroPosteriorImbalance(obj)
      forceSC = obj.SC.forceVector;
      forceIS = obj.IS.forceVector;
      forceTM = obj.TM.forceVector;
      scapulaCS = obj.shoulder.scapula.coordSys;

      forceResultant = Vector(scapulaCS.origin, scapulaCS.origin) + ...
        forceSC + forceIS + forceTM; 

      % Take infero-superior orthogonal component <=> project on scapular axial plane
      inferoSuperior = Vector(scapulaCS.origin, scapulaCS.origin + scapulaCS.IS);
      imbalanceVector = inferoSuperior.orthogonalComplementToPoint(forceResultant.target);

      lateroMedial = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.ML);
      anteroPosterior = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.PA);
      obj.anteroPosteriorImbalance = ...
        sign(dot(imbalanceVector, anteroPosterior)) *...
        rad2deg(lateroMedial.angle(imbalanceVector));
    end
  end

end