function measureUniquePointInsertionEquivalent(obj, humerus)
  % Find and return the point for which the sidemost muscle's segmentation 
  % centroids are coplanar with muscle's insertion endpoints and humerus center.
  % This unique point insertion equivalent can then be used with existing
  % methods to compute the muscle's fibres contact points. 

  if size(obj.centroid, 1) < 2
    obj.uniquePointInsertionEquivalent =...
      getUniqueCentroidInsertionEquivalent(obj, humerus);
    return 
  end

  insertionEndpoints =...
    humerus.landmarks.insertions.endpoints.(obj.getFullName())';
  insertionDirection =...
    Vector(insertionEndpoints(1,:), insertionEndpoints(2,:));

  sidemostCentroids = getSidemostCentroids(obj.centroid, insertionDirection);

  sidemostPlanesNormal = [];
  for i = 1:size(sidemostCentroids, 1)
    sidemostPlane = Plane;
    sidemostPlane.fit(...
      [humerus.center; sidemostCentroids(i,:); insertionEndpoints(i,:)]);
    sidemostPlanesNormal = [sidemostPlanesNormal; sidemostPlane.normal];
  end

  sidemostPlanesIntersection =...
    Vector(cross(sidemostPlanesNormal(1,:), sidemostPlanesNormal(2,:)));

  humerusCenterToInsertionEndpoints =...
    Vector(humerus.center, mean(insertionEndpoints));

  humerusCenterToInsertionEquivalent =...
    sidemostPlanesIntersection.orientToward(humerusCenterToInsertionEndpoints);

  insertionEquivalent =...
    humerus.center...
    + humerus.radius * humerusCenterToInsertionEquivalent.direction;

  obj.uniquePointInsertionEquivalent = insertionEquivalent;
end

function output = getUniqueCentroidInsertionEquivalent(muscle, humerus)
  % In this case the unique point insertion equivalent is chosen to be in the
  % middle of insertion endpoints
  humerusCenterToInsertionEquivalent = Vector(humerus.center,...
    mean(humerus.landmarks.insertions.endpoints.(muscle.getFullName()))' );

  output = humerus.center...
    + humerus.radius * humerusCenterToInsertionEquivalent.direction;
end

function output = getSidemostCentroids(centroids, insertionDirection)
  % This algorithm is sufficient to determine the sidemost centroids given that
  % all the centroids are coplanar
  centroidsPositionAlongInsertion = [];
  for i = 1:size(centroids, 1)
    centroidsPositionAlongInsertion(i) = insertionDirection.dot(centroids(i,:));
  end

  sidemostCentroids = [];
  sidemostCentroids = [sidemostCentroids; centroids(find(...
    centroidsPositionAlongInsertion == min(centroidsPositionAlongInsertion)),:)];
  sidemostCentroids = [sidemostCentroids; centroids(find(...
    centroidsPositionAlongInsertion == max(centroidsPositionAlongInsertion)),:)];

  output = sidemostCentroids;
end