clc
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                                                       %
%   Skript soužící pro:     extrakci tlouštky z Metadat programu OctsegCopy,                                            %    
%                       
%                                                                                                                       %
%   Nutne soubory:  ID_scan_number.vol - OCT volume scan                                                                %
%                   ID_scan_number.JPG - fundus snímek                                                                  %                      
%                                                                                                                       % 
%   ve vsech metadatech provedeny manulani korekce zajmovych oblasti.                                                   %                       
%   Potreba zvolit:    adresar - cesta bez posledniho podadresare                                                       % 
%                      ID_scan_number - nazev slozky, kde jsou ulozeny Meta data a zaroven nazev metadat                %     
%                      
%                                                                                                                       %
%   Pokud není vytvořena tlouťková mapa, je spuštěna extrakce této mapy. Při tom je potřeba rucne zvolit oblast         %
%   optickeho disku. V teto oblasti je nastavna nulová tloušťka NFL. Skript vytvori soubor s nazvem:                    %  
%   thickness_ID_scan_number.mat.                                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all; clc;

%% Volitelne položky:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
adresar = 'C:\Users\Fierka\Documents\BAKALARKA\data\normalni\'; % cesta k datum bez posledniho podadresare 
ID_scan_number='71930' % nazev slozky, kde jsou ulozeny Meta data a zaroven musi byt s stejny korenovy nazev metadat 
addpath ('C:\Users\Fierka\Documents\BAKALARKA\data\normalni\'); % Cesta do podadresáře s programem pro registraci
okno=20; % Volba velikosti okna pro výběr vzorku z fundus snímku (vel. 20 = okno 41x41 px)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

okno_OCT = round(okno/2.35); %  cislo na konci slouzi ke korekci rozdilnosti hustoty pixelu u fundus a SLO (slo by nejak automatizovat???)
addpath ([adresar, ID_scan_number '\']); % Cesty k datům OCT+fundus
DatabasePath = [adresar, ID_scan_number '\'];

%% Načtení souborů:
OCTname = [ID_scan_number '.vol']; % Název OCT volume scanu             
FUNname = [ID_scan_number '.JPG']; % Název fundus snímku

%% OCT data loading

% fid = fopen(['c:\DISERTACE\Fundus\programy\markus_OCT\octsegCopy\data\' ID_scan_number '\' OCTname]);
fid = fopen([adresar, ID_scan_number '\' ID_scan_number '.vol' ]);

% File header read
header.Version = fread( fid, 1, 'int8' );
header.SizeX = fread( fid, 1, 'int32' );
header.NumBScans = fread( fid, 1, 'int32' );
header.SizeZ = fread( fid, 1, 'int32' );
header.ScaleX = fread( fid, 1, 'double' );
header.Distance = fread( fid, 1, 'double' );
header.ScaleZ = fread( fid, 1, 'double' );
header.SizeXSlo = fread( fid, 1, 'int32' );
header.SizeYSlo = fread( fid, 1, 'int32' );
header.ScaleXSlo = fread( fid, 1, 'double' );
header.ScaleYSlo = fread( fid, 1, 'double' );
header.FieldSizeSlo = fread( fid, 1, 'int32' );
header.ScanFocus = fread( fid, 1, 'double' );
header.ScanPosition = fread( fid, 4, 'schar' );
header.ExamTime = fread( fid, 2, 'int32' );
header.ScanPattern = fread( fid, 1, 'int32' );
header.BScanHdrSize = fread( fid, 1, 'int32' );
header.ID = fread( fid, 16, 'int8' );
header.ReferenceID = fread( fid, 16, 'int8' );
header.PID = fread( fid, 1, 'int32' );
header.PatientID = fread( fid, 21, 'int8' );
header.Padding = fread( fid, 3, 'int8' );
header.DOB = fread( fid, 1, 'double' );
header.VID = fread( fid, 1, 'int32' );
header.VisitID = fread( fid, 24, 'int8' );
header.VisitDate = fread( fid, 1, 'double' );
header.Spare = fread( fid, 1840, 'int8' );

disp(['           Version: ' char(header.Version')]); 
disp(['             SizeX: ' num2str(header.SizeX)]);
disp(['         NumBScans: ' num2str(header.NumBScans)]);
disp(['             SizeZ: ' num2str(header.SizeZ)]);
disp(['            ScaleX: ' num2str(header.ScaleX) ' mm']);
disp(['          Distance: ' num2str(header.Distance) ' mm']);
disp(['            ScaleZ: ' num2str(header.ScaleZ) ' mm']);
disp(['          SizeXSlo: ' num2str(header.SizeXSlo)]);
disp(['          SizeYSlo: ' num2str(header.SizeYSlo)]);
disp(['         ScaleXSlo: ' num2str(header.ScaleXSlo) ' mm']);
disp(['         ScaleYSlo: ' num2str(header.ScaleYSlo) ' mm']);
disp(['FieldSizeSlo (FOV): ' num2str(header.FieldSizeSlo) '°']);
disp(['         ScanFocus: ' num2str(header.ScanFocus) ' dpt']);
disp(['      ScanPosition: ' char(header.ScanPosition')]); %OD - right eye; OS - left eye
% disp(['          ExamTime: ' datestr(header.ExamTime(1)/(1e7*60*60*24)+1601*365.25)]);
disp(['       ScanPattern: ' num2str(header.ScanPattern)]);
disp(['      BScanHdrSize: ' num2str(header.BScanHdrSize) ' bytes']);
disp(['                ID: ' char(header.ID)']);
disp(['       ReferenceID: ' char(header.ReferenceID')]);
disp(['               PID: ' num2str(header.PID)]);
disp(['         PatientID: ' char(header.PatientID')]);
disp(['               DOB: ' datestr(header.DOB+693975)]);
disp(['               VID: ' num2str(header.VID)]);
disp(['           VisitID: ' char(header.VisitID')]);
disp(['         VisitDate: ' datestr(header.VisitDate+693975)]);


% status = fseek( fid, 2048, -1 );


%% SLO image read
slo = fread( fid, header.SizeXSlo*header.SizeYSlo, 'uint8' );
slo = reshape( slo, header.SizeXSlo, header.SizeYSlo );
slo = imrotate( slo, 90 );
slo = flipud( slo );


%% B-scans read (including image data and scanning position)
% status = fseek( fid, 2048+(header.SizeXSlo*header.SizeYSlo), -1 );

BScans=zeros(header.SizeZ,header.SizeX,header.NumBScans);
ILM=zeros(header.NumBScans-1,header.SizeX);
RPE=zeros(header.NumBScans-1,header.SizeX);

if header.NumBScans == 1
    NScans=1;
else
    NScans=header.NumBScans-1;  
end

%Evaluation of BScan position in SLO image
ii = 1;
        status = fseek( fid, 16+2048+(header.SizeXSlo*header.SizeYSlo)+(ii*(header.BScanHdrSize+header.SizeX*header.SizeZ*4)), -1 );
        StartX1 = fread( fid, 1, 'double' );
        StartX1_px=round(StartX1/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %StartX1 in pixels
        StartY1 = fread( fid, 1, 'double' );
        StartY1_px=round(StartY1/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %StartY1 in pixels
        EndX1 = fread( fid, 1, 'double' );
        EndX1_px=round(EndX1/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %EndX1 in pixels
        EndY1 = fread( fid, 1, 'double' );
        EndY1_px=round(EndY1/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %EndY1 in pixels
%         NumSeg = fread( fid, 1, 'int32' );
        
ii = NScans;
        status = fseek( fid, 16+2048+(header.SizeXSlo*header.SizeYSlo)+(ii*(header.BScanHdrSize+header.SizeX*header.SizeZ*4)), -1 );
        StartX2 = fread( fid, 1, 'double' );
        StartX2_px=round(StartX2/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %StartX2 in pixels
        StartY2 = fread( fid, 1, 'double' );
        StartY2_px=round(StartY2/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %StartY2 in pixels
        EndX2 = fread( fid, 1, 'double' );
        EndX2_px=round(EndX2/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %EndX2 in pixels
        EndY2 = fread( fid, 1, 'double' );
        EndY2_px=round(EndY2/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %EndY2 in pixels
        NumSeg = fread( fid, 1, 'int32' );

% clear line

%% Fundus data loading

fun=double(imread(FUNname));
R=fun(:,:,1);
G=fun(:,:,2);
B=fun(:,:,3);
funGB=(G+B)/2;


%% Zjisteni, zda uz byl vytvoren soubor s tloustkovou mapou a soubor s registraci fundus+OCT

DatabaseFiles = dir(DatabasePath);

ValidThicknessFile = [];
ValidRegistrationFile = [];

for i = 1:size(DatabaseFiles,1)
    if regexpi(DatabaseFiles(i).name, 'thickness_');
        ValidThicknessFile = 1;
    end
    if regexpi(DatabaseFiles(i).name, 'reg_param_');
        ValidRegistrationFile = 1;
    end
end


if ValidThicknessFile
    load(['thickness_' ID_scan_number '.mat']); 

else 
%% Vypocet tloustky z Meta souborů 
NumberOfScanes = 2;
    f_name= ['thickness_' ID_scan_number '.mat'];

    for i = 1:size(DatabaseFiles,1)
        ValidMetaFile = regexpi(DatabaseFiles(i).name, '.meta');
        if ValidMetaFile
            NumberOfScanes = NumberOfScanes + 1; % Number of all meta files
            mode='DefaultMedlineAutoData'; % nutno upravit
            typemode = 'num';
            output_IN = [];
            output_OUT = [];
            fid = fopen([DatabasePath,'\',DatabaseFiles(i).name], 'r');
            if fid == -1
                return;
            end

            linie = fgetl(fid);

            while (numel(output_IN) == 0) && ischar(linie)
                [des rem] = strtok(linie);
                if strcmp(des, mode)
                    if strcmp(typemode, 'str')          
                        output_IN = rem;
                    else
                        output_IN = str2num(rem);
                    end
                    break; 
                end
                linie = fgetl(fid);
            end
            
            mode = 'DefaultRPEautoData'; % nutno upravit

            while (numel(output_OUT) == 0) && ischar(linie)
                [des rem] = strtok(linie);
                if strcmp(des, mode)
                    if strcmp(typemode, 'str')          
                        output_OUT = rem;
                    else
                        output_OUT = str2num(rem);
                    end
                    break; 
                end
                linie = fgetl(fid);
            end
        %%%%%%
        NFL_thick = output_OUT-output_IN; 
        NFL_thick(NFL_thick<0)=0; % Pokud je INFL výš než ONFL, zvol nulovou tloušťku 
        NFL_thick(NFL_thick>50)=50; % Pokud je tloušťka větší než 50, jde patrně o chybu....
        NFL_thickness_RAW(:,NumberOfScanes) = NFL_thick;
        output_IN_RAW(NumberOfScanes,:)=output_IN
        output_OUT_RAW(NumberOfScanes,:)=output_OUT 
        
        %%%%%%
         max_IN(:,NumberOfScanes)=max(max(output_IN_RAW(NumberOfScanes,:))); % pro OPL_RPE
         min_OUT(:,NumberOfScanes)=min(min(output_OUT_RAW(NumberOfScanes,:))); % pro OPL_RPE

        end
        ValidMetaFile = [];
    end
    fclose(fid);
    
    if NumberOfScanes == 0
        disp('Nejsou ulozeny MetaData, je potreba nejprve spustit soubor octsegMain.m a oznacit vrstvy')
        return
    end
    
    NFL_thickness_RAW = NFL_thickness_RAW';
    NFL_thickness_RAW = flipdim(NFL_thickness_RAW,1);
    NFL_thickness_RAW = round(NFL_thickness_RAW);

    
    EndY2_px(EndY2_px<1)=1;
    StartY1_px(StartY1_px<1)=1;
    StartX1_px(StartX1_px<1)=1;
    EndX1_px(EndX1_px<1)=1;
    SLO_resized = slo(EndY2_px:StartY1_px,StartX1_px:EndX1_px); 

    [mm,nn] = size(SLO_resized);

    NFL_thickness_resized = imresize(NFL_thickness_RAW, [mm nn]);


   save ([DatabasePath 'thickness_' ID_scan_number '.mat'],'max_IN','min_OUT','output_OUT_RAW','output_IN_RAW','NFL_thickness_RAW', 'NFL_thickness_resized'); % pro OPL_RPE
    disp(' ')
    disp('Byl vytvorel tloustkovy soubor')
end
