%function [fileSet] = loadFile(filename)
function [fileSet, sortedFile] = loadFile(filename)
% function loadFile enables the user to load a file with patient data 
%and sorts the slices according to the DICOM information 
%into separate uniform volume sets

%the output: 
%>sortedFile . . . is a cell structure containing names of the slices
% acquired from DICOMinfo information and sorted according to their
% z-coordinate (also present); additional information about spacing and size 
% is included (this output can be omitted)

%>fileSet. . . is a cell containing the processed volumes; first
% column is implicitly set to be the result according to the defined
% function and the second column represents the result as a parametric image
% of one eigenvalue. The third column stores the voxel size information and
% the next is the origin of the volume

%loading of the non-single z-slice sequences exclusively
filename = strcat(filename, '/');
imName =strcat(filename, '00000000.dcm');
k = 0;j = 1;
while exist(imName)
    infoSlice = dicominfo(imName);
    if (isfield(infoSlice,'ImagesInAcquisition')~=0) && ((abs(infoSlice.ImagePositionPatient(3)) - abs(infoSlice.SliceLocation))<0.01)
            slices{j,3} = infoSlice.SliceLocation;
            slices{j,1} = imName;
            slices{j,2} = [(infoSlice.PixelSpacing)' infoSlice.SpacingBetweenSlices];
        j = j+1;
    end 
    k = k+1;
    imNumber =strcat('0000', num2str(k));
    while length(imNumber)<8
        imNumber = strcat('0',imNumber);
    end        
    imName = strcat(filename,imNumber,'.dcm');
end

% Gaussian mask setting for preprocessing
Mask = 3;                                                                                               
G = fspecial('gaussian', [Mask Mask],1);

%recursive function that sorts the loaded sequences with additional
%simplification of the sequence storage
sortedFile= volumes(slices,1,{});
sS = size(sortedFile);
for i=1:sS(2)
    volSorted =  sortedFile{1,i};
    spacing = volSorted{1,2};
    volSorted(:,2) = [];
    sortedFile{1,i} = volSorted;
end

%inicialization of the ouput cell
fileSet = [];
for i=1:sS(2)
    vol = sortedFile{1,i};
    sV = size(vol);
    spacing(i,:) = spacing;
    origin = vol{1,2};
    orgVol(:,:,1)=dicomread(vol{1,1});
    NoSame = 1;k=1;
    % the images are loaded into a 3Dmatrix, if they are at the same position, they are averaged
    for j=2:sV(1)  
        j
        if (vol{j,2})==(vol{j-1,2})
            NoSame = NoSame+1;
            orgVol(:,:,k) = (double(dicomread(vol{j,1}))+((NoSame-1)*(double(orgVol(:,:,k)))))/NoSame;           
        else k=k+1
             orgVol(:,:,k)=dicomread(vol{j,1});
             NoSame=1;
        end
    end
    [l,w,v] = size(orgVol);
    %smoothening with the Gaussian mask
    for j=1:v
    orgVolF(:,:,j) = imfilter(orgVol(:,:,j),G, 'same');
    end   
 vol=[];orgVol=[]; 

% inicialization of the output matrices 

 lambda = zeros(l,w,v);
 result = zeros(l,w,v);

%  Preparation of the convolution masks for the computation of the Hessian
%  matrix
for k =1:1
    h{k}=cat(k,[1],[-2],[1]);    
  spaceMat(k) = 1/(spacing(i,k)^2);
end
h{4}=[1 0 -1;0 0 0; -1 0 1];
h{5}=cat(3,[1 0 -1]',[0 0 0]',[-1 0 1]');
h{6}=cat(3,[1 0 -1],[0 0 0],[-1 0 1]);       
 spaceMat(4) = 1/(spacing(i,1)*spacing(i,2));
 spaceMat(5) = 1/(spacing(i,1)*spacing(i,3));
 spaceMat(6) = 1/(spacing(i,2)*spacing(i,3));
 spaceMat(4:6) = 1/4 * spaceMat(4:6);

    for j=1:6
    mask = spaceMat(j) * h{j};
    H{j}=convn(orgVolF,mask,'same');
    end
%  The Hessian matrix values
    H11=H{1}; H22=H{2}; H33=H{3};
    H12=H{4};H13=H{5}; H23=H{6};
    
%    The computation of the eigenvalues
    for j=1:(l*w*v)
    mat = [H11(j) H12(j) H13(j);H12(j) H22(j) H23(j);H13(j) H23(j) H33(j)];
    val = eig(mat);
    sortVal=sort(val);
    lambda(j) =abs(sortVal(1));
    l1(j) =sortVal(2);
    l2(j) =sortVal(3);
    result(j) =sqrt(l1(j).^2+l2(j).^2);
    end 
%Assigning the result and normalizing the volume, also storing values as unsigned char

resultF = result-min(min(min(result)));
resultF = resultF/max(max(max(resultF)));
resultF = round(resultF*255); 

lambdaF = lambda-min(min(min(lambda)));
lambdaF = lambdaF/max(max(max(lambdaF)));
lambdaF = round(lambdaF*255);    
% Storing the processed volume into a cell structure
fileSet{i,2}=resultF;
fileSet{i,3}=lambdaF;
fileSet{i,1}=spacing(i,:);
fileSet{i,4}= origin;
orgVolF=[];lambdaF=[];resultF=[];
end


end





