close all; clear all;clc;format longg;
warning off;
%########################################################################
%
%	match_lenstar_to_ultrasound
%
%	This code was written in 2018 by Mustapha El Hamdaoui and Rafael Grytz
%
% 	This code comes with ABSOLUTELY NO WARRANTY. This is free
%	software, and you are welcome to redistribute it.
%
%	The code was tested using MATLAB version R2016b and R2018a.
%
%########################################################################
%
% Please modify the following variables when using the code with your data:
%
% Specify the sound the sound velocities of [v_CT v_ACD v_LT v_VCD]
velocity = [1.610 1.540 1.7233 1.540];
%
% Specify the center wavelength of your LenStar in nm
lambda = 831.78;
%
% Specify the percentage of data intended for the fitting.
% We used 75% for training (fitting) and 25% for validation.
train = 75;
%
% Specify if you want to use a random subset (rand_data=true) or the 
% top portion (rand_data=false) of the data for the fitting.
% To reproduce the results in our paper use rand_data=false.
rand_data = false;
%
% In this code, we assume that CT and ACD can (cannot) be resolved form 
% the ultrasound measures if the number of columns of the data file is 
% even (odd).
%
%########################################################################
%
% tol is the fitting tollerance
tol   = 1e-8;
% n0 is the initial guess used to minimize the SSR (fitting)
n0    = 1.5;
% maximum number of iterations
max_iter = 100;
% nCT is the human corneal refractive index estimated at the center
% wavelength lambda.
nCT = 1.361594+3*6009.687/lambda^2-5*676076000/lambda^4+7*59084500000000/lambda^6;
% After you run the script, a window will pop up and allow you to
% navigate to the file that contains the data
FileNAme     = uigetfile('*.xls*');
% Reading the chosen file
 data     = xlsread(FileNAme);
%##########################################################################

l_data = length(data);
if train < 100
    l_train=floor(train*length(data)/100);
    if rand_data
        % Use random portion of data for fitting
        rnd = randperm(l_data);
        data_train      = data(rnd(1:l_train),:);
        data_val        = data(rnd(l_train+1:end),:);
    else
        % Use top portion of data for fitting
        data_train   = data(1:l_train,:);
        data_val     = data(l_train+1:end,:);
    end
    % We save both data_train and data_valid as ascii files
    save -ascii data_fitting.txt data_train
    save -ascii data_validation.txt data_val
else
    data_train = data;
end
%
% Sorting data intp OPLs and GLs
if mod(size(data,2),2) == 0
    comp_num = size(data,2)/2;
    iplots   = comp_num;
    iplots0  = 0;
    for i  = 1:comp_num
        OPL(:,i) = data_train(:,i);
        GL(:,i)  = data_train(:,i+comp_num) * velocity(i);
        OPL_val(:,i) = data_val(:,i);
        GL_val(:,i)  = data_val(:,i+comp_num) * velocity(i);
    end
    comps = {'CT' 'ACD' 'LT' 'VCD'};
else
    comp_num = (size(data,2)-3)/2;
    iplots   = comp_num+1;
    iplots0  = 1;
    for i  = 1:comp_num
        OPL(:,i) = data_train(:,i+2);
        GL(:,i)  = data_train(:,i+3+comp_num) * velocity(i+2);
        OPL_val(:,i) = data_val(:,i+2);
        GL_val(:,i)  = data_val(:,i+3+comp_num) * velocity(i+2);
    end
    t_ASD(:)       =  data_train(:,3+comp_num);
    GL_CT(:)       =  data_train(:,1)/nCT;
    OPL_ACD(:)     =  data_train(:,2);
    t_ASD_val(:)   =  data_val(:,3+comp_num);
    GL_CT_val(:)   =  data_val(:,1)/nCT;
    OPL_ACD_val(:) =  data_val(:,2);
    comps = {'LT' 'VCD'};
end

% Fitting matching data sets
ft              = fittype('x/n','independent','x','dependent','y');
opts            = fitoptions('Method','NonlinearLeastSquares');
opts.Robust     = 'Bisquare';
opts.StartPoint = n0;
opts.TolFun     = tol;
opts.TolX       = tol;
RI              = zeros(1,size(OPL,2));
R2              = zeros(1,size(OPL,2));
d               = 0:0.01:10;
for i = 1:size(OPL,2)
    [xData,yData]  = prepareCurveData(OPL(:,i),GL(:,i));
    [result,gof]   = fit(xData,yData,ft,opts);
    RI(i)          = coeffvalues(result);
    R2(i)          = gof.rsquare*100;
    F_m            = d/RI(i);
    h = figure(1);
    subplot(1,iplots,iplots0+i)
    plot(xData,yData,'*r', d,F_m,'k')
    xlabel(sprintf('OPL_{LS,%s} [mm]',comps{i}));
    ylabel(sprintf('GL_{US,%s} [mm]',comps{i}));
    legend('Data Ponts','Best Fit','Location','southeast')
    axis([min(xData)-min(xData)/20 max(xData)+max(xData)/20 min(yData)-min(yData)/20 max(yData)+max(yData)/20])
    title(strcat('RI=', num2str(RI(i)),' ; R^2=',num2str(ceil(R2(i))),'%'))
end

if mod(size(data,2),2) ~= 0
    err   = 10;
    iter  = 0;
    ratio_0 = 0.8;
    while (abs(err) > tol) && (iter < max_iter)
        vASD      =  velocity(1) / (1 + (ratio_0/velocity(2))*(velocity(1)-velocity(2)));
        GL_ASD(:) = t_ASD(:) * vASD;
        GL_ACD    = GL_ASD-GL_CT;
        [xData,yData] = prepareCurveData(OPL_ACD,GL_ACD);
        [result, gof] = fit( xData, yData, ft, opts );
        nACD    = coeffvalues(result);
        Rs      = gof.rsquare*100;
        ratio   = mean(OPL_ACD./(nACD*GL_ASD));
        err     = (ratio - ratio_0)/ ratio;
        ratio_0 = ratio;
        iter    = iter+1;
    end
    LT_est      = data_val(:,3)/RI(1);
    F_m_ACD = d/nACD;
    subplot(1,iplots,1)
    plot(OPL_ACD,GL_ACD,'*r',d,F_m_ACD,'k')
    xlabel('OPL_{LS,ACD} [mm]');
    ylabel('GL_{ACD} [mm]');
    legend('Data Ponts','Best Fit','Location','southeast')
    axis([min(xData)-min(xData)/20 max(xData)+max(xData)/20 min(yData)-min(yData)/20 max(yData)+max(yData)/20])
    title(strcat('RI=', num2str(nACD),' ; Rs=',num2str(ceil(Rs)),'%'))
    fprintf('v_ASD = %.6f\n',vASD);
    fprintf('GL_ACD/GL_ASD = %.6f\n',ratio);
    fprintf('n_ACD = %.6f\n',nACD);
end
for i = 1:size(OPL,2)
    fprintf('n_%s = %.6f\n',comps{i},RI(i));
end
set(h,'Position',[0 50 300*(iplots+1) 300])

% Generate box plot using validation data set
if train < 100
    b_plot = figure(2);
    AL_LS = 0.0;
    AL_US = 0.0;
    for i = 1:size(OPL_val,2)
        box_err(:,i) =  100*(OPL_val(:,i)/RI(i) - GL_val(:,i))./GL_val(:,i);
        AL_LS = AL_LS + OPL_val(:,i)/RI(i);
        AL_US = AL_US + GL_val(:,i);
    end
    
    if mod(size(data,2),2) ~= 0
        ASD_LS = GL_CT_val(:) + OPL_ACD_val(:)/nACD;
        ASD_US = t_ASD_val(:) * vASD;
        AL_LS = AL_LS + ASD_LS;
        AL_US = AL_US + ASD_US;
        ASD_err = 100*(ASD_LS-ASD_US)./ASD_US;
        box_err = [ASD_err box_err];
        comps = ['ASD', comps];
    end
    AL_err = 100*(AL_LS-AL_US)./AL_US;
    box_err = [box_err AL_err];

    boxplot(box_err,'Labels',[comps, 'AL'])
    ylabel('Error [%]')
    grid on
    set(b_plot,'Position',[0 450 600 500]);
end

