%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%David B. Berry || dbberry@ucsd.edu                                       %
%This script will calculate:                                              %
%      -Angle with Respect to the Horizontal                              %
%      -Sagittal Cobb Angle                                               %
%      -Coronal Cobb Angle                                                %
%      -Posterior IVD Height                                              %
%      -Central IVD Height                                                %
%      -Anterior IVD Height                                               %
%      -Sagittal IVD Angle                                                %
%      -Coronal IVD Angle                                                 %
%      -Degree of Axial Rotation                                          %
%      -Upper Lumbar Sagittal Cobb Angle                                  %
%      -Lower Lumbar Sagittal Cobb Angle                                  %
%                                                                         %
%                                                                         %
%The inputs for the script are:                                           %
%      -.csv ROIs on "corners" of vertebrae in sagittal images            %
%      -.csv ROIs on the posterior elements of axially reformatted imgs   %
%                                                                         %
%                                                                         %
%The this script calls the following functions which must be in the       %
% /Functions Directory                                                    %
%      -filenamer                                                         %
%      -readROI2                                                          %
%      -redROI2postvec                                                    %
%      -FitPlaneDB                                                        %
%      -coeffcheck                                                        %
%      -rotatecoords                                                      %
%      -getAngle                                                          %
%      -derotate                                                          %
%      -getHorizontalangleDB                                              %
%      -getIVDdistances                                                   %
%      -vectorplotter2                                                    %
%      -titleheader                                                       %     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all; clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Please change these parameters   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ipres=1.2; %in plane resolution mm^2
thk=5; %slice thickness mm
num_subjects=1;
num_positions=1;
displayflag = 1; %Display images representing the ROIs & planes?
writedataflag = 1; %Write data to a .txt file?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  First locate .csv files of ROIs  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[~, path]=uigetfile('*.csv',...
    'Please select a Sagittal file in the directory you want to analyze');
[~, path2]=uigetfile('*.csv',...
    'Please select an Axial file in the directory you want to analyze');

tic %Start time
addpath('Functions')
addpath(path)
addpath(path2)
files = dir([path '/*.csv']);
files_Axial = dir([path2 '/*.csv']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Error message to ensure same       %
%number of CSV files in each folder %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(files,1) ~= size(files_Axial,1)
    error('FileSizeError',...
        'There is not an equal number of sagittal and axial .csv''s')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Begin loop to calculate parameters.%
%Each iteration calculates one      %
%position for each subject.         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(files)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Extract ROI coordinates and scale %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf('Analyzing file: %0.f\n', i)
    [file, file_ax]=filenamer(files,i);%Get Axial/Sagittal filenames
    [x_coord,y_coord,z_coord, ~, ~, ~]=...
        readROI2(file);%Read and sort Sagittal ROIs
    [postvec]=readROI2postvec(file_ax,ipres,thk);%Calculate the x vector
    %from the ROIs that defines the posterior elements
    x_coord=x_coord*ipres; y_coord=y_coord*thk; z_coord=z_coord*ipres;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %NOTE: To this point I have stripped%
    %the data from the csv file. It is  %
    %organized to x_coord & x_coord_ped %
    %(expand for y, z).  The prior gives%
    %all of the  coordinates for the n  %
    %dimension for the vertebral        %
    %endplates, the later defines the   %
    %coordinates for the pedicle base.  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Fit planes to superior and inferior%
    %endplates                          %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [snormal,scoeff,splanes,sCenter]=...
        FitPlaneDB(x_coord,y_coord,z_coord,'S',0);
    [inormal,icoeff,iplanes,iCenter]=...
        FitPlaneDB(x_coord,y_coord,z_coord,'I',0);
    syaxis=cell2mat(snormal)';
    iyaxis=cell2mat(inormal)';
    scoeff=coeffcheck(scoeff);
    icoeff=coeffcheck(icoeff);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %NOTE: Load x axis & the component  %
    %of the x axis || to it from the    %
    %coeff matrix (2nd) column.  Then   %
    %figure out the rotation between the%
    %two of them.                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Assign the Z and X axes            %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    [rotsyaxis]=rotatecoords(snormal, postvec);
    [rotiyaxis]=rotatecoords(inormal, postvec);
    for ll=1:length(rotiyaxis)
        if inormal{ll}(3) < -0
            rotiyaxis{ll}=[-rotiyaxis{ll}(1,:);rotiyaxis{ll}(2,:);...
                -rotiyaxis{ll}(3,:)];
        end
        if snormal{ll}(3) < -0
            rotsyaxis{ll}=[-rotsyaxis{ll}(1,:);rotsyaxis{ll}(2,:);...
                -rotsyaxis{ll}(3,:)];
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Switch from Osirix defs to ISB defs%
    %Osirix: Y Right Z Caudal X Anterior%
    %ISB: Y Caudal Z Left X Anterior    %
    %X=nxaxis; Z=rotyaxis; Y=zaxis;     %
    %NOTE:20131024 Z in the matrix is Y %
    %Y is Z, X is X.  Everything is     %
    %transposed Y=nxaxis; Zs=rotsyaxis; %
    %Xs=zsaxis;                         %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %       Calculate Euler Angles      %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for j=1:length(rotsyaxis)-1
        [psi(j), theta(j), phi(j)]=getAngle(rotiyaxis{j},rotsyaxis{j+1});
        [psi(j), theta(j), phi(j)]=...
            derotate(rotiyaxis{j},rotsyaxis{j+1},psi(j), theta(j), phi(j));
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate Distance Between Origins %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [x_trans,y_trans,z_trans]=...
        origin_translation(x_coord, y_coord, z_coord);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   Calculate Sagittal Cobb Angle   %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [cpsi, ctheta, cphi]=getAngle(rotsyaxis{1},rotsyaxis{6});
    [cpsi, ctheta, cphi]=...
        derotate(rotsyaxis{1},rotsyaxis{6},cpsi, ctheta, cphi);
    sag_cobb=[ctheta zeros(1,4)];
    cor_cobb=[cpsi zeros(1,4)];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %        Get Horizontal Angle       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    horiz=[getHorizontalAngleDB(x_coord,z_coord) zeros(1,4)];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %       Calculate IVD Heights       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for w=1:length(splanes)-1
        [DAnterior(w), DCenter(w), DPosterior(w)] = getIVDdistances(0,0,...
            splanes{w+1}(:,1:5),splanes{w+1}(:,6:10),...
            splanes{w+1}(:,11:15),iplanes{w}(:,1:5),iplanes{w}(:,6:10),...
            iplanes{w}(:,11:15));
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate Upper Sagittal Cobb Angle%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [cpsiu, cthetau, cphiu]=getAngle(rotsyaxis{1},rotiyaxis{3});
    [cpsiu, cthetau, cphiu]=...
        derotate(rotsyaxis{1},rotiyaxis{3},cpsiu, cthetau, cphiu);
    cobbu=[cthetau zeros(1,4)];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate Lower Sagittal Cobb Angle%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [cpsil, cthetal, cphil]=getAngle(rotsyaxis{4},rotsyaxis{6});
    [cpsil, cthetal, cphil]=...
        derotate(rotsyaxis{4},rotiyaxis{5},cpsil, cthetal, cphil);
    cobbl=[cthetal zeros(1,4)];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %         Figures to display        %
    %Vectorplotter gives sup and inf    %
    %coordinates superimposed for each  %
    %level. FitPlaneDB gives plane & ROI%
    %for each level in a 3D plot.       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if displayflag
        vectorplotter2(rotsyaxis,rotiyaxis)
        figure('Color','w')
        [~,~,~]=FitPlaneDB(x_coord,y_coord,z_coord,'S',displayflag);
        [~,~,~]=FitPlaneDB(x_coord,y_coord,z_coord,'I',displayflag);
        axis equal
        displayflag=1; %Plot only once when analyzing many files
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Store data in cell array to write %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fw=files(i).name(1:5);
    pp=['P' num2str(i-1)];
    subject={fw fw fw fw fw}';
    position={pp pp pp pp pp}';
    level={'L1-L2' 'L2-L3' 'L3-L4' 'L4-L5' 'L5-S1'}';
    A=strread(num2str(theta),'%s');
    B=strread(num2str(-1*psi),'%s');
    C=strread(num2str(phi),'%s');
    D=strread(num2str(DAnterior),'%s');
    E=strread(num2str(DCenter),'%s');
    F=strread(num2str(DPosterior),'%s');
    G=strread(num2str(sag_cobb),'%s');
    H=strread(num2str(cor_cobb),'%s');
    I=strread(num2str(cobbu),'%s');
    J=strread(num2str(cobbl),'%s');
    K=strread(num2str(horiz),'%s');
    writemat{i,:}=[subject position level A B C D E F G H I J K];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  Clear variables for next sample  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear x_coord y_coord z_coord x_coord_ped y_coord_ped z_coord_ped
    clear postvec
    clc
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Write data to .txt file      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if writedataflag==1
    xlmat={'Subject'; 'Position'; 'Level'; 'Theta (Cor)'; 'Psi (Sag)';...
        'Phi (Axi)'; 'DAnterior'; 'DCenter'; 'DPosterior';...
        'Sagittal Cobb Angle';'Coronal Cobb Angle';...
        'Upper Cobb Angle';'Lower Cobb Angle';'Horizontal Angle'}';
    for e=1:length(writemat)
        xlmat=[xlmat; writemat{e}];
    end

    xlmat2=titleheader(num_positions);
    
    for J=1:num_subjects
        matwrite{J,1}=xlmat{20*J-20+2,1};
        for I=1:num_positions
            matwrite{J,2+35*I-35}=xlmat{20*J-20+2+I*5-5,10};
            matwrite{J,3+35*I-35}=xlmat{20*J-20+2+I*5-5,11};
            matwrite{J,4+35*I-35}=xlmat{20*J-20+2+I*5-5,12};
            matwrite{J,5+35*I-35}=xlmat{20*J-20+2+I*5-5,13};
            matwrite{J,6+35*I-35}=xlmat{20*J-20+2+I*5-5,14};
            for K=1:5
                matwrite{J,7+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,4};
                matwrite{J,8+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,5};
                matwrite{J,9+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,6};
                matwrite{J,10+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,7};
                matwrite{J,11+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,8};
                matwrite{J,12+I*35-35+(K*6-6)}=xlmat{20*J-20+1+I*5-5+K,9};
            end
        end
    end
    xlmat=[xlmat2; matwrite];
    ds=cell2dataset(xlmat);
    export(ds,'file',['Spine_measurements' date '.txt']) %Write .txt file
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Can be opened with excel by right  %
%clicking the file and select       %
%open\with->Excel.                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc, toc %display total script time