function [GC, varargout] = glenoid_orientation_preparation(implant_type, x_s, y_s, z_s, SN, GlenoidImplant_Center, image_position, varargin)

% -------------------------------------------------------------------------
% Description
% -------------------------------------------------------------------------
% This function defines the glenoid fossa center, glenoid version, and
% glenoid inlcination in the scapula coordinate system to be used in the
% subject-specific tool of the MSM in order to replicat patinets with
% anatomical total shoudelr arthroplasty.
% The values are defined according to
% "Terrier, A., et al. "Measurements of three-dimensional glenoid erosion
% when planning the prosthetic replacement of osteoarthritic shoulders.",
% The bone & joint journal 96.4 (2014): 513-518." and also the description
% provided in documaent named "MSM_input_output_preprocessing".
%
% -------------------------------------------------------------------------
% Syntax:
% Case 1: for anatomical total shoulder arthroplasty
% [GC, GV, GI] = glenoid_orientation_preparation('aTSA', x_s, y_s, z_s, SN, GlenoidImplant_Center, image_position, GlenoidImplant_LateralAxis)
% GC = glenoid_orientation_preparation(rTSA, x_s, y_s, z_s, SN, GlenoidImplant_Center, image_position)

% -------------------------------------------------------------------------
% by Ehsan Sarshari, summer 2018 
% -------------------------------------------------------------------------
% clear, clc
% 
% % set the current path
% current_folder = pwd;
% 
% % differentiate between platforms
% if isunix
%     data_file = [current_folder ,'/', 'patients_data.mat'];
% else
%     data_file = [current_folder ,'\', 'patients_data.mat'];
% end
% 
% % load the patients data
% load(data_file)
% 
% % define the number of patinets for who you uploaded the data
% number_of_patients = length(fieldnames(patients_data));
% 
% % ids of the patients
% patients_ids = fieldnames(patients_data);
% 
% % initialaize the results
% GV_S = zeros(number_of_patients,1);
% GI_S = zeros(number_of_patients,1);
% GC_S = zeros(number_of_patients,3);

%     
% for patient_counter = 1: number_of_patients
%    
%     
%     % load the scapula coordinate axes
%     x = eval(['patients_data.', patients_ids{patient_counter}, '.X_S']);
%     y = eval(['patients_data.', patients_ids{patient_counter}, '.Y_S']);
%     z = eval(['patients_data.', patients_ids{patient_counter}, '.Z_S']);


    % define the rotation matrix from scpula coordinate system to CT coordinate
    % system
    R_S_CT = [x_s', y_s', z_s'];
    
    if strcmp(implant_type, 'aTSA') % if it's natomical
        
        GlenoidImplant_LateralAxis = varargin{:};

        % define the GC_axis in the scapula frame
        GC_axis_S = R_S_CT'*GlenoidImplant_LateralAxis';

        % define the version and inlcination in the scapula coordinate system
        % (consistent with the MSM definition)
        GV = atand(GC_axis_S(1)/GC_axis_S(3)); varargout{2} = GV;% glenoid version

        GI = atand(GC_axis_S(2)/GC_axis_S(3)); varargout{3} = GI;% glenoid inlcination
        
    else
        % it's rTSA and we do not have GI or GV for it.
%         GV = 'NA';
%         GI = 'NA';
    end
    

    % define the GC in CT by adding the image position to the GC from
    % blueprint data
    GC_CT = GlenoidImplant_Center' + image_position;
    
    % define the glenoid fossa center in the scapula coordinate system
    GC = (R_S_CT'*(GC_CT - SN') + [0 0 4]')';    
    
return

% % save the data
% if isunix
%     save_file = [current_folder ,'/', 'MSM_inputs.mat'];
% else
%     save_file = [current_folder ,'\', 'MSM_inputs.mat'];
% end
% 
% save(save_file, 'GV_S', 'GI_S', 'GC_S');