clc
close all
clear all

%% Nacteni tloustkoveho souboru 
adresar = 'C:\Users\Fierka\Documents\BAKALARKA\data\normalni\'; % cesta k datum bez posledniho podadresare 
ID_scan_number='71084' % 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 podadresare s programem pro registraci
addpath ([adresar, ID_scan_number '\']); % Cesty k datum OCT+fundus
DatabasePath = [adresar, ID_scan_number '\'];

DatabaseFiles = dir(DatabasePath);

ValidThicknessFile = [];

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


if ValidThicknessFile
    load(['thickness_' ID_scan_number '.mat']); 
else 
    disp ('Neni vytvoren tloustkovy soubor!!!')
end
%% Soubor pro ktery se tvori 2D projekce
path=(['C:\Users\Fierka\Documents\BAKALARKA\data\normalni\' ID_scan_number '\' ID_scan_number '.vol']); %uprava pro skript
show=1; 


fid = fopen(path);

%% Nacteni hlavicky souboru
header.Version = fread( fid, 12, '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 );


%% Nacteni SLO obrazu
slo = fread( fid, header.SizeXSlo*header.SizeYSlo, 'uint8' );
slo = reshape( slo, header.SizeXSlo, header.SizeYSlo );
slo = imrotate( slo, 90 );
slo=slo/max(max(slo)); % double

% figure(1)
% imshow(slo,[])
% imwrite(slo,'slo.png')
% slo = flipud( slo );
% figure(1); subplot(1,2,1);
% imshow(slo,[]); 
% hold on
% scrsz = get(0,'ScreenSize');
% figure('Position',[1 0 scrsz(3) scrsz(4)-70]);
% clf;
% if show==1
%     figure(2);subplot(1,2,1);
%     imshow(slo,[]);
% end

%% Nacteni B-skenu (zahrnuje obrazova data a pozici skenovani) 
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

    for ii = 1:NScans
%         ii
        % Hodnoceni pozice B-skenu v SLO obraze
        status = fseek( fid, 16+2048+(header.SizeXSlo*header.SizeYSlo)+(ii*(header.BScanHdrSize+header.SizeX*header.SizeZ*4)), -1 );
        StartX = fread( fid, 1, 'double' );
        StartX_px=round(StartX/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %PocatekX v pixelech
        StartY = fread( fid, 1, 'double' );
        StartY_px=round(StartY/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %PocatekY in pixelech
        EndX = fread( fid, 1, 'double' );
        EndX_px=round(EndX/(header.SizeXSlo*header.ScaleXSlo)*header.SizeXSlo); %KonecX in pixelech
        EndY = fread( fid, 1, 'double' );
        EndY_px=round(EndY/(header.SizeYSlo*header.ScaleYSlo)*header.SizeYSlo); %KonecY in pixelech
        NumSeg = fread( fid, 1, 'int32' );
%         subplot(1,2,1);
%         if show==1
%             figure(2);
%             subplot(1,2,1); line([StartX_px EndX_px],[StartY_px EndY_px]);
%         end
        
        %Nacteni B-skenu
        status = fseek( fid, header.BScanHdrSize+2048+(header.SizeXSlo*header.SizeYSlo)+(ii*(header.BScanHdrSize+header.SizeX*header.SizeZ*4)), -1 );
        oct = fread( fid, header.SizeX*header.SizeZ, 'float32' );
        oct = power( reshape( oct, header.SizeX, header.SizeZ ), 0.25 );
        oct = flipud( oct );
        oct = imrotate( oct, -90 );
        BScans(:,:,ii)=oct;
       
%         subplot(1,2,2);
%         if show ==1
%             figure(2);
%             subplot(1,2,2); imshow( BScans(:,:,ii) );
%             drawnow
%         end
    % Nacitani segmentovanych dat
        status = fseek( fid, 256+2048+(header.SizeXSlo*header.SizeYSlo)+(ii*(header.BScanHdrSize+header.SizeX*header.SizeZ*4)), -1 );
        seg = (fread( fid, NumSeg*header.SizeX, 'float' ))';
        ILM(ii,:) = seg(1:header.SizeX);
        RPE(ii,:) = seg(header.SizeX+1:header.SizeX*2);
%         pause 

radky(ii)=[StartY_px]; %vektor obsahujici pozice na nichz se v SLO obrazu nachazi jednotlive B-skeny
    end
    fclose( fid );
%% Vlastni komentovany skript pro vytvoreni 2D projekci
%% Vypocteni prumeru
output_OUT=round(output_OUT_RAW);
AA=size(output_IN_RAW);
for k=1:header.NumBScans 
    oct_sken=BScans(:,:,k); % vezme dany B-sken
for l=3:AA(1)
    for m=1:header.SizeX   
    prumer(k,m)=mean(oct_sken(output_IN_RAW(l,m):output_OUT(l,m),m)); 
        % pro dany B-sken vypocte prumer z hodnot mezi pozicemi ohranicenymi pri segmentaci
    prumery=[prumer];
    end
end
end

%% Vytvoreni 2D projekce v rozmerech 
projekce2Dprumer=zeros((radky(1)-radky(end)), header.SizeX); %vytvoreni matice o velikosti snimku
for o=1:NScans
    projekce2Dprumer((radky(o)-(radky(end)-1)),:)=[prumery(o,:)]; %doplneni hodnot na pozice B-skenu
end

S=size(projekce2Dprumer);
for q=1:header.SizeX 
    for p=1:(S(1)-1) 
if projekce2Dprumer(p,q)~=0 %prohledava matici dokud nenalezne nenulovou hodnotu
x1=p;               % radek, na kterem se nachazi nenulova hodnota
y1=projekce2Dprumer(p,q); % nenulova hodnota
index=p+1;          % pozice sousedici s nenulovou hodnotou
while projekce2Dprumer(index, q)==0&&index<S(1)
                    % dokud jsou nulove, prochazi hodnoty na sousednich pozicich
                    % a zaroven je tato pozice stale v matici
index=index+1;      
end
x2=index;           % radek, na kterem se nachazi dalsi nenulova hodnota
y2=projekce2Dprumer(index,q); % nenulova hodnota na dane pozici
a=interp1([x1, x2], [y1, y2], linspace( x1, x2, (x2-x1+1)), 'spline'); % vypocet interpolaci mezi nenulovymi hodnotami
projekce2Dprumer(x1:x2,q)=a; % interpolované hodnoty se doplni na 'prazdna' mista v matici
end
    end
end
imtool(projekce2Dprumer)

%% Vypocteni Medianu
output_OUT=round(output_OUT_RAW);
AA=size(output_IN_RAW);
for k=1:header.NumBScans 
    oct_sken=BScans(:,:,k); % vezme dany B-sken
for l=3:AA(1)
    for m=1:header.SizeX   
    med(k,m)=median(oct_sken(output_IN_RAW(l,m):output_OUT(l,m),m));
            % pro dany B-sken vypocte median z hodnot mezi pozicemi
            % ohranicenymi pri segmentaci
    meds=[med];
    end
end
end

%% Vytvoreni 2D projekce v rozmerech 
projekce2Dmedian=zeros((radky(1)-radky(end)), header.SizeX); %vytvoreni matice o velikosti snimku
for o=1:NScans
    projekce2Dmedian((radky(o)-(radky(end)-1)),:)=[meds(o,:)]; %doplneni hodnot na pozice B-skenu
end

S=size(projekce2Dmedian);
for q=1:header.SizeX 
    for p=1:(S(1)-1) 
if projekce2Dmedian(p,q)~=0 %prohledava matici dokud nenalezne nenulovou hodnotu
x1=p;               % radek, na kterem se nachazi nenulova hodnota
y1=projekce2Dmedian(p,q); % nenulova hodnota
index=p+1;          % pozice sousedici s nenulovou hodnotou
while projekce2Dmedian(index, q)==0&&index<S(1)
                    % dokud jsou nulove, prochazi hodnoty na sousednich pozicich
                    % a zaroven je tato pozice stale v matici
index=index+1;      
end
x2=index;           % radek, na kterem se nachazi dalsi nenulova hodnota
y2=projekce2Dmedian(index,q); % nenulova hodnota na dane pozici
a=interp1([x1, x2], [y1, y2], linspace( x1, x2, (x2-x1+1)), 'spline'); % vypocet interpolaci mezi nenulovymi hodnotami
projekce2Dmedian(x1:x2,q)=a; % interpolované hodnoty se doplni na 'prazdna' mista v matici
end
    end
end
imtool(projekce2Dmedian)
    %% Nalezeni Minima
output_OUT=round(output_OUT_RAW);
AA=size(output_IN_RAW);
for k=1:header.NumBScans 
    oct_sken=BScans(:,:,k); % vezme dany B-sken
for l=3:AA(1)
    for m=1:header.SizeX   
    minim(k,m)=min(oct_sken(output_IN_RAW(l,m):output_OUT(l,m),m)); 
        % pro dany B-sken vypocte prumer z hodnot mezi pozicemi ohranicenymi pri segmentaci
    minima=[minim];
    end
end
end
%% Vytvoreni 2D projekce v rozmerech 
projekce2Dminima=zeros((radky(1)-radky(end)), header.SizeX); %vytvoreni matice o velikosti snimku
for o=1:NScans
    projekce2Dminima((radky(o)-(radky(end)-1)),:)=[minima(o,:)]; %doplneni hodnot na pozice B-skenu
end

S=size(projekce2Dminima);
for q=1:header.SizeX 
    for p=1:(S(1)-1) 
if projekce2Dminima(p,q)~=0 %prohledava matici dokud nenalezne nenulovou hodnotu
x1=p;               % radek, na kterem se nachazi nenulova hodnota
y1=projekce2Dminima(p,q); % nenulova hodnota
index=p+1;          % pozice sousedici s nenulovou hodnotou
while projekce2Dminima(index, q)==0&&index<S(1)
                    % dokud jsou nulove, prochazi hodnoty na sousednich pozicich
                    % a zaroven je tato pozice stale v matici
index=index+1;      
end
x2=index;           % radek, na kterem se nachazi dalsi nenulova hodnota
y2=projekce2Dminima(index,q); % nenulova hodnota na dane pozici
a=interp1([x1, x2], [y1, y2], linspace( x1, x2, (x2-x1+1)), 'spline'); % vypocet interpolaci mezi nenulovymi hodnotami
projekce2Dminima(x1:x2,q)=a; % interpolované hodnoty se doplni na 'prazdna' mista v matici
end
    end
end
imtool(projekce2Dminima); %zobrazeni 2D projekce
%% Nalezeni Maxima
output_OUT=round(output_OUT_RAW);
AA=size(output_OUT);
for k=1:header.NumBScans 
    oct_sken=BScans(:,:,k); % vezme dany B-sken
for l=3:AA(1)
    for m=1:header.SizeX   
    maxim(k,m)=max(oct_sken(output_IN_RAW(l,m):output_OUT(l,m),m));
     % pro dany B-sken vezme maximum z hodnot mezi pozicemi ohranicenymi pri segmentaci
    maxims=[maxim];
    end
end
end
%% Vytvoreni 2D projekce v rozmerech 
projekce2Dmaxima=zeros((radky(1)-radky(end)), header.SizeX); %vytvoreni matice o velikosti snimku
for o=1:NScans
    projekce2Dmaxima((radky(o)-(radky(end)-1)),:)=[maxims(o,:)]; %doplneni hodnot na pozice B-skenu
end

S=size(projekce2Dmaxima);
for q=1:header.SizeX 
    for p=1:(S(1)-1) 
if projekce2Dmaxima(p,q)~=0 %prohledava matici dokud nenalezne nenulovou hodnotu
x1=p;               % radek, na kterem se nachazi nenulova hodnota
y1=projekce2Dmaxima(p,q); % nenulova hodnota
index=p+1;          % pozice sousedici s nenulovou hodnotou
while projekce2Dmaxima(index, q)==0&&index<S(1)
                    % dokud jsou nulove, prochazi hodnoty na sousednich pozicich
                    % a zaroven je tato pozice stale v matici
index=index+1;      
end
x2=index;           % radek, na kterem se nachazi dalsi nenulova hodnota
y2=projekce2Dmaxima(index,q); % nenulova hodnota na dane pozici
a=interp1([x1, x2], [y1, y2], linspace( x1, x2, (x2-x1+1)), 'spline'); % vypocet interpolaci mezi nenulovymi hodnotami
projekce2Dmaxima(x1:x2,q)=a; % interpolované hodnoty se doplni na 'prazdna' mista v matici
end
    end
end
imtool(projekce2Dmaxima)
disp('Vsechnz projekce byly vytvoreny')
