function [I,Iorg]=SortDicomFiles_OUH(DirPath,FileNameFilter,recursive,pathToDCMTK,AvoidTerminalOutput)
  %This function scans a directory (resurively if resursive=true) and sort
  %the located dicom files acording to Study, Series and Modality
  %The returned value Iorg is identical to the result that is obtained from
  %the function FindDicomFileNames. The result I has a structure similar to
  %Iorg however where Iorg has the filename of the dicom file I has the entire
  %dicom header. The function is almost the same speed as FindDicomFiles
  %but the returned Dicomheaders are usefull in a number of situations and
  %can potentially be used to speed up the reading of the files
  %Please see specification of FileNameFilter in the function RecursiveDir_OUH
    %19-03-2015 CaB
  
  %Program will use the function ReadLargePinnacleStructFile_OUH to read
  %struct file in case the struct file is too large for matlab to read (bug in relation to matlab)
  %2015-09-15 CAB
  
  %Change the program to allways read RTSTRUCT files of a size larger than
  %1 Mebabyte by the program ReadLargePinnacleStructFile_OUH. This new
  %version need access to the dcmdump program from dmctk
  %2015-09-22 CAB
  
  %Change back to version 2015-0915. The problem was not the code
  %but execution of the code on Matlab2014 machines whithout the needed DICOM patch
  %2015-12-08 CAB
  
  % Added option to give path to dicom toolkit as parameter, to be used by
  % ReadLargePinnacleStructFile_OUH. Removed check for dicom toolkit in
  % this function, as well as old versions of code commented out to improve
  % readability. Old versions are found in SVN.
  % 2016-02-26 RST
  
  %11-07-2016 changed dicominfo call to use 'UseVRHeuristic',false which is
  %a recomended feature if the redaing the data fails. This change seems to
  %correct the problem reading pinnacle struct files. /CaB
  
  %Have made a change such that RTPLANs and RTSTRUCT are not
  %grouped together even if their study and series UID match (e.g. normal
  %and escelated plan based on the same CT dataset)
  %31-08-2016 CaB
  
  %Have made a change such that MR images with different diffusion B-values
  %are not grouped together even if their study and series UID match.
  %19-01-2017 UBe
  
  %Have changes the spatial sorting of the files such that it also work for
  %sagital images
  %24-02-2017 CaB
  
  %Have changed such that dose files with same seris and study uid are not
  %grouped togehter if they are poiting to different plans. A problem seen
  %in Oncentra
  %25-09-2017 CaB
  %Smalle correction to the above that it is validated that if
  %ImageOrientationPAtient exsit it should have a length of 6 (non-empty)
  %03-03-2019 CaB
  
  %dicominfo has caused older versions of Matlab to crash when reading 
  %large structure sets from Pinnacle. A new function, GetModality_OUH, 
  %is introduced to use dcmtk to determine which files need to be 
  %processed with special care.
  %2019-12-03 SLK
  
  %An initial scan for doublicated UID has been added. Only one file per
  %UID will be handled in the code (could be introduced by duplicated data
  %in the data directory)
  %2020-03-24 CaB
  
  %Added input parameter (AvoidTerminalOutput), which, if set to true, turns off output to the terminal
  %Also added that an input of an empty string for pathToDCMTK is equivalent to effectively not providing that input
  %The default values for the individual parameters are:
  %FileNameFilter='.'
  %recursive=true
  %pathToDCMTK=""
  %AvoidTerminalOutput=false
  %2023-01-16 CaB

  %Corrected a bug that occurred for CT images, which does not have the tag Manufacture defined.
  %In that case, the previous version did not detect the images, which is updated with the current bug fix.
  %2024-08-28 CaB

  if nargin<1
    return;
  end
  if nargin<2
    FileNameFilter='.'; %all files are included
  end
  if nargin<3
    recursive=true;
  end
  if nargin>=4
    if pathToDCMTK==""
      clear('pathToDCMTK')
    end
  end
    
  if nargin<5
    AvoidTerminalOutput=false;
  end
  files=RecursiveDir_OUH(DirPath,FileNameFilter,recursive);
  Data=cell(size(files,1),1);
  counter=1;
  for i=1:size(files,1)
    if ~AvoidTerminalOutput
      if mod(i,ceil(size(files,1)/100))==0
        temp=i/size(files,1);
        disp(temp);
      end
    end
    try
      %From version 2015-09-15 
      lastwarn('');
      if verLessThan('matlab', '8.6')
        if strcmp(GetModality_OUH(files{i}, pathToDCMTK),'RTSTRUCT')
            Data{counter}=ReadLargePinnacleStructFile_OUH(files{i},pathToDCMTK);
        else
            Data{counter}=dicominfo(files{i});
        end
      else
        Data{counter}=dicominfo(files{i},'UseVRHeuristic', false);
      end
      if ~isempty(strfind(lastwarn,'Not enough data imported.'))
         lastwarn('');
         if ~exist('pathToDCMTK','var')
            Data{counter}=ReadLargePinnacleStructFile_OUH(files{i});
         else
            Data{counter}=ReadLargePinnacleStructFile_OUH(files{i},pathToDCMTK);
         end
      end
      %Since we are only interested in DicomData with attached information
      %(e.g. images struct plan dose) where are some dicom files which we
      %might disregard e.g. DicomDirectory files. Below is there are possibility
      %to list the MediaStorageSOPClassUID (SOP class) which we do not like to load.
      %The value 1.2.840.10008.1.3.10 is related to DicomDirectory files
      %and 1.3.46.670589.11.0.0.12.1; 1.3.46.670589.11.0.0.12.2; 1.3.46.670589.11.0.0.12.4 are Philips private SOP classes
      %2017-07-26 CAB
      AcceptData=true;
      try
        temp=Data{counter}.MediaStorageSOPClassUID;
        SOPNotToLoad={'1.2.840.10008.1.3.10','1.3.46.670589.11.0.0.12.1','1.3.46.670589.11.0.0.12.2','1.3.46.670589.11.0.0.12.4'};
        if sum(cellfun(@(x) strcmp(x,temp),SOPNotToLoad))>0.5
          AcceptData=false;
        end
      catch
      end
      if AcceptData
        counter=counter+1;
      end
    catch e
      warning(['The file: ' files{i} ' is not a DICOM file \n' e.message]) 
    end
  end
  Data=Data(1:counter-1);
  
  %Have added a scan for UID before sorting all the Dicom images. If the
  %same UID is detected, only one of these files will be evaluated. This
  %scan is made to handle the situation for which the same Dicom file is
  %present more than once within the given directory. Previously the scan
  %of UID has performed within each Grouped set of images. That method had
  %the consequence that two plans or struct files would be grouped in
  %different groups, and even if they had the same UID, the data would be
  %handled for both identical files.
  %CaB 2020-03-24
  UIDS=cell(length(Data),1);
  for i=1:length(Data)
    try
      UIDS{i}=Data{i}.SOPInstanceUID;
    catch
      UIDS{i}=[];
      warning('Some of the files do not have a UID value');
    end
  end
  [~,ia,~]=unique(UIDS);
  ia=unique([ia ; find(cellfun(@(x) isempty(x),UIDS))]);
  Data=Data(ia);
  clear UIDS;
    
  FileNamesDicomInfo=cell(9,2);
  FileNamesDicomInfo{1,1}='RTPLAN';   FileNamesDicomInfo{1,2}=cell(0,0); 
  FileNamesDicomInfo{2,1}='RTDOSE';   FileNamesDicomInfo{2,2}=cell(0,0);
  FileNamesDicomInfo{3,1}='CT';       FileNamesDicomInfo{3,2}=cell(0,0);
  FileNamesDicomInfo{4,1}='RTSTRUCT'; FileNamesDicomInfo{4,2}=cell(0,0);
  FileNamesDicomInfo{5,1}='MR';       FileNamesDicomInfo{5,2}=cell(0,0);
  FileNamesDicomInfo{6,1}='PT';      FileNamesDicomInfo{6,2}=cell(0,0);
  FileNamesDicomInfo{7,1}='NM';       FileNamesDicomInfo{7,2}=cell(0,0);
  FileNamesDicomInfo{8,1}='CBCT';     FileNamesDicomInfo{8,2}=cell(0,0);
  FileNamesDicomInfo{9,1}='NA';       FileNamesDicomInfo{9,2}=cell(0,0);
  
  
  if ~isempty(Data)
  %Check existence of Study and Series UID and Modality
  index=cellfun(@(x) ~isfield(x,'StudyInstanceUID'),Data);
  if sum(index)~=0
    warning('Some of the files does not contain a Study UID. Will be grouped together');
    temp=1:size(index,1);
    for i=temp(index)
      Data{i}.StudyInstanceUID='';
    end
  end
  index=cellfun(@(x) ~isfield(x,'SeriesInstanceUID'),Data);
  if sum(index)~=0
    warning('Some of the files does not contain a Series UID. Will be grouped together');
    temp=1:size(index,1);
    for i=temp(index)
      Data{i}.SeriesInstanceUID='';
    end
  end
  index=cellfun(@(x) ~isfield(x,'Modality'),Data);
  if sum(index)~=0
    warning('Some of the files does not contain Modality info. Will be grouped as NA');
    temp=1:size(index,1);
    for i=temp(index)
      Data{i}.Modality='NA';
    end
  end
  Study_Series=cellfun(@(x) [x.StudyInstanceUID '_' x.SeriesInstanceUID],Data,'UniformOutput',false);
  
  %Avoid to group RTplan and RTStruct togehter even if there study and series
  %UID are the same
  index=cellfun(@(x) strcmpi(x.Modality,'RTPLAN')||strcmpi(x.Modality,'RTSTRUCT') ,Data);
  for i=find(index==1)'
    Study_Series{i}=[Study_Series{i} dicomuid];
  end
  
  %Avoid to group MR images with different diffusion B-values even if the
  %study and series UID are the same
  index=cellfun(@(x) isfield(x,'DiffusionBValue'),Data);
  for i=find(index==1)'
    Study_Series{i}=[Study_Series{i} '_' num2str(Data{i,1}.DiffusionBValue)];
  end
  %Avoid to group images which have the main component of the vector orthogonal
  %to images plane in different x y z direction even if the study and series UID
  %are the same. This avoid grouping Surveys with different normal vector into
  %the same group
  %2017-07-26 CAB
  %Have also added that if ImageOritationPatient exist is also need to be
  %of length 6 (thus non-empty)
  %2019-03-03 CAB
  index=cellfun(@(x) isfield(x,'ImageOrientationPatient') && length(x.ImageOrientationPatient)==6,Data);
  for i=find(index==1)'
    temp=cross(Data{i,1}.ImageOrientationPatient(1:3),Data{i,1}.ImageOrientationPatient(4:6));
    temp=find(abs(temp)==max(abs(temp)));
    temp=temp(1);
    Study_Series{i}=[Study_Series{i} '_' num2str(temp)];
  end
  %Avoid to group dose together even if they do have the same series and
  %study UID but is referencing different plan UID. This has been observed
  %for some versions of Oncentra
  %2017-09-25 CAB
  index=cellfun(@(x) strcmpi(x.Modality,'RTDOSE'),Data);
  for i=find(index==1)'
    %Study_Series{i}=[Study_Series{i} dicomuid];
    if isfield(Data{i},'ReferencedRTPlanSequence')
      if isfield(Data{i}.ReferencedRTPlanSequence,'Item_1')
        if isfield(Data{i}.ReferencedRTPlanSequence.Item_1,'ReferencedSOPInstanceUID')
          temp=Data{i}.ReferencedRTPlanSequence.Item_1.ReferencedSOPInstanceUID;
          Study_Series{i}=[Study_Series{i} '_' temp];
        end
      end
    end
  end
  
  ListOfStudySeries=unique(Study_Series);
  DataGrouped=cell(size(ListOfStudySeries,1),1);
  for i=1:size(ListOfStudySeries,1)
    index=cellfun(@(x) strcmp(x,ListOfStudySeries{i}),Study_Series);
    DataGrouped{i}=Data(index);
  end

  for i=1:size(DataGrouped,1)
    indextot=false(size(DataGrouped{i},1),1);
    for j=1:size(FileNamesDicomInfo,1)-2
      index=cellfun(@(x) strcmp(x.Modality,FileNamesDicomInfo{j,1}),DataGrouped{i});
      indextot=indextot | index;
      if sum(index)~=0
        if strcmp(FileNamesDicomInfo{j,1},'CT')
          %Check if it is a CT or CBCT scan
          if isfield(DataGrouped{i}{1},'Manufacturer')
            if strncmpi(DataGrouped{i}{1}.Manufacturer,'ELEKTA',6)||strncmpi(DataGrouped{i}{1}.Manufacturer,'Varian',6) 
              FileNamesDicomInfo{8,2}{end+1}=DataGrouped{i}(index);
            else
              FileNamesDicomInfo{j,2}{end+1}=DataGrouped{i}(index);  
            end
          else
            FileNamesDicomInfo{j,2}{end+1}=DataGrouped{i}(index);
          end
        else
          FileNamesDicomInfo{j,2}{end+1}=DataGrouped{i}(index);
        end
      end
    end
    if sum(~indextot)~=0
      FileNamesDicomInfo{9,2}{end+1}=DataGrouped{i}(~indextot);
    end
  end
  
  %The next section has been removed at the same time that the initial scan
  %for duplicated UID was added earlier in the code.
  %2020-03-24 CaB
  
  %Remove images with identical UID (e.g. same dataset exist twice in the
  %folder)
%   for j=1:size(FileNamesDicomInfo,1)
%     if ~isempty(FileNamesDicomInfo{j,2})
%       for i=1:size(FileNamesDicomInfo{j,2},2)
%         try
%           UIDS=cellfun(@(x) x.SOPInstanceUID,FileNamesDicomInfo{j,2}{i},'UniformOutput',false);
%         catch
%           UIDS=[];
%         end
%         if ~isempty(UIDS)
%           [~,ia,~]=unique(UIDS);
%           FileNamesDicomInfo{j,2}{i}=FileNamesDicomInfo{j,2}{i}(ia);
%         else
%           warning('Some of the files do not have a UID value');
%         end
%       end
%     end
%   end
  %Sort CT/CBCT/MR/PET images according to position
  for j=[3 5 6 8]
    for i=1:size(FileNamesDicomInfo{j,2},2)
      try
        %Find first index for position sorting. Normally the third index
        %but not if sagital images are used. temp is a absolute vector orthorgonal
        %to the image plane except for the use of abs
        temp=cellfun(@(x) cross(x.ImageOrientationPatient(1:3),x.ImageOrientationPatient(4:6)),FileNamesDicomInfo{j,2}{i},'UniformOutput',false);
        temp=cellfun(@(x) find(abs(x)==max(abs(x))),temp,'UniformOutput',false);
        temp=cellfun(@(x) x(1),temp);
        indexImagePos=temp(1);
        temptest=temp==indexImagePos;
        if ~all(temptest)
          warning('Direction orthogonal to image plan seem to change for the images');
        end
        positions=cellfun(@(x) x.ImagePositionPatient(indexImagePos),FileNamesDicomInfo{j,2}{i});
      catch
        positions=[];
      end
      if ~isempty(positions)
        [~,ix]=sort(positions);
        FileNamesDicomInfo{j,2}{i}=FileNamesDicomInfo{j,2}{i}(ix);
      end
    end
  end
  %Time sort all entries
  for j=1:size(FileNamesDicomInfo,1)
    if ~isempty(FileNamesDicomInfo{j,2})
      DateNumber=zeros(size(FileNamesDicomInfo{j,2},2),1);
      for i=1:size(FileNamesDicomInfo{j,2},2)
        try
          Time=FileNamesDicomInfo{j,2}{i}{1}.StudyTime;
          if ~ischar(Time)  || length(Time)~=6
            Time='000001';
          end
          Date=FileNamesDicomInfo{j,2}{i}{1}.StudyDate;
          if ~ischar(Date)  || length(Date)~=8
            Date='19000101';
          end
        catch
          Time='000001';
          Date='19000101';
        end
        
        DateNumber(i) = datenum([Date Time],'yyyymmddHHMMSS');
       
      end
      [~,ix]=sort(DateNumber);
      FileNamesDicomInfo{j,2}=FileNamesDicomInfo{j,2}(ix);
    end
  end
  I=FileNamesDicomInfo;
  Iorg=I;
  for i=1:size(Iorg,1)
    for j=1:size(Iorg{i,2},2)
      Iorg{i,2}{j}=cellfun(@(x) x.Filename ,Iorg{i,2}{j},'UniformOutput',false);
    end
  end
  else
    I=FileNamesDicomInfo;
    Iorg=I;
  end
end

