%   This function computes the Back Vertex Power and Addition of a 
%   Multizone Contact Lens (MCLs) by means of the exported Power Profile
%   (PP) from NIMO TR1504
%
%   REQUIRED RESOURCES 
%
%   (1) Perform a PP measure with NIMO TR1504 and export the 
%   "Averaged Power" in a .CSV file.
%   (2) This function uses the vline function developed by Brandon Kuczenski
%   which can be downloaded from Matlab File Exchange:
%   https://es.mathworks.com/matlabcentral/fileexchange/1039-hline-and-vline
%   If this function is not downloaded and added to the path the
%   nimoProcess function will return an error when the parameter graph is
%   set to 1.
%
%   PARAMETERS
%
%   zones       Integer number which indicates the zone pairs 
%               that are going to be considered in order to compute the BVP 
%               and Add (i.e. for MCLs with 6 alternating near and distance zones
%               the total number of 3 zone pairs should be inserted. 
%
%   ignoreCZ    Ignore Central Zone: Optional double number which indicates the central 
%               zone (radial position from the center in mm) of the lens that is 
%               not going to be considered during the zone recognition proccess. 
%               This parameter is required in MCLs that presents a central defect which
%               may be confussed with a theoretical zone. 
%
%   sens        Double number which indicates the sensitivity of 
%               the recognition zone process in order to detect what changes 
%               in power (D) correspond with the real transition of a zone. 
%               Lower values detect lower slopes in the zones transition.
%               The recommended value by the authors is 0.2 D.
%
%   graph       Binary number which indicates if a plot with the power
%               profile is going to be presented with vertical lines
%               indicating the zones transition. (0 for not presenting
%               graph and 1 for presenting).
%   
%   RETURN VARIABLES 
%
%   BVP         Back Vertex Power of the MCLs in Diopters (D)
%   Add         Addition of the lens
%   
%   function [BVP,Add]= nimoProcess(zones,ignoreCZ,sens,graph)
%
%   EXAMPLE
%   
%   Process a PP of a MCLs of 6 alternating zones considering that there is
%   a central defect that creates a false zone in the central 0.3 mm of
%   diameter considering minimal changes of power of 0.2D among zones. Show
%   the graph with the power profile.
%
%   function [BVP,Add]= nimoProcess(3,0.3,0.2,1)

function [BVP,Add]=nimoProcess(zones,ignoreCZ,sens, graph)

%------- IMPORTING DATA FROM CSV FILE TO MATLAB ---------%

%Open dialog box for selecting the exported *.csv file from NIMO.
[FileName,PathName]=uigetfile({'*.csv'},'Select the exported CSV file from NIMO'); 

%"data" variable is a two columns matrix which contains the radial position
%in the first column at the corresponding power (D) in the second column.
data=csvread(strcat(PathName,'/',FileName),1,0);

%------- DETECT ZONES LOCATION ---------%

x=data(:,1);
y=data(:,2);
%Computing the first derivate
first_derivate=y(1:end-1,1)-y(2:end); 
%Inflection points
a=abs(abs(first_derivate)>sens); 
%Computing the second derivate
second_derivate=a(1:end-1,1)-a(2:end); 
%Select the limit points
point=find(second_derivate,length(second_derivate),'first'); 
%Location of the limit points that separate the zones.
Order=point+1; 
limit_points=x(Order);
%Draws the graph
if graph
figure,plot (x,y),vline(x(Order))
end
%Detect the optical axis of the MCL (center).
zeroposition=length(data)/2+1;
%Detect the vector position corresponding to the ignoreCZ
[minim,cutposition]=min(abs(data(zeroposition:end,1)-ignoreCZ)); 
%Find the nearest value in the data to omit the central zone.
prox=data(cutposition+zeroposition); 
%The locations of the inflection points in the data for the positive radial
%positions.
limit_points=[cutposition+zeroposition;Order(length(Order)/2+1:end)]; 

%------- PROCESSING OF THE INFORMATION OF THE ZONES ---------%

%Zone sizes for al the recognized zones including transition area.
size_zones=limit_points(2:end)-limit_points(1:end-1);
%M is a matrix of NaN defined by size_zones variable 
M=NaN(max(size_zones),length(size_zones));
%Inserting power values for each one of the zones. Each columns represents
%a zone starting by the first zone in the first column.
for i=1:length(size_zones)
M(1:size_zones(i),i)=data(limit_points(i):(limit_points(i+1)-1),2);
end

%------- COMPUTING BVP AND ADD VALUES ---------%
%Reshape of the matrix M in a vector that contains the values of the Power
Pb=reshape(M(:,1:4:length(size_zones)),1,prod(size(M(:,1:4:length(size_zones)))))'; 
%Average of the baseline power values.
Zone1=nanmean(Pb(1:max(size_zones)*zones)); 
%Reshape of the matrix M in a vector that contains the values of the
%addition.
Padition=reshape(M(:,3:4:length(size_zones)),1,prod(size(M(:,3:4:length(size_zones)))))';
%Average of all points of the addition zone.
Zone2=nanmean(Padition(1:max(size_zones)*zones));
%Computing the addition power
Add=Zone2-Zone1; 
%Computing the Back Vertex Power
BVP=Zone1; 

disp('The value of the Addition is: ');
disp(Add);
disp('The value of the Back Vertex Power is: ');
disp(BVP);
end