%%
% This function sorts all DICOM files from a given directory tree into 
% seperate directories based on study and series information. Input is the 
% root of the directory tree containing the files. Output is the root of the 
% directory tree under which the files are copied to.
%
% The script provides similar functionality to tht provided by the DicomSorter
% program written by Jolinda Smith (Lewis Center for NeuroImaging, 
% University of Oregon). The reason for reimplementing the same 
% functionality is that I only have a windows executable of that
% program. This script is "multi platform" - runs on all platforms MATLAB
% does. Also the availability of source code allows others to modify according
% to their needs.
%
% REQUIREMENT: Image processing toolbox.
%
% WARNING: Function will overwrite data if there are name conflicts. 
%
% may 2014: bug fix for a bug reported by He Su - if output directory is a
%           sub-directory of the input directory then the script will
%           process the input multiple times.
%
% @author: Ziv Yaniv (zyaniv@childrensnational.org)
%

%
% Copyright Ziv Yaniv
% 
%    Licensed under the Apache License, Version 2.0 (the "License");
%    you may not use this file except in compliance with the License.
%    You may obtain a copy of the License at
% 
%      http://www.apache.org/licenses/LICENSE-2.0
% 
%    Unless required by applicable law or agreed to in writing, software
%    distributed under the License is distributed on an "AS IS" BASIS,
%    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%    See the License for the specific language governing permissions and
%    limitations under the License.
%
function dicomSorter(inputDirectoryName, outputDirectoryName)

        %check that we received the path to two existing directories
if ~(isdir(inputDirectoryName) && isdir(outputDirectoryName))
  error('Invalid input: input or output directory does not exist.');
end

        %check that the ouput directory is not a subdirectory of the input
        %this will lead to duplicates after some files are already copied
        %to the output directory (when we reach it during traversal of the
        %directory structure these "output" files become "input")
fullInputPath = cd(cd(inputDirectoryName));
fullOutputPath = cd(cd(outputDirectoryName));
if strfind(fullOutputPath, fullInputPath)
  error('Output directory cannot be a subdirectory of the input directory.');
end
       %this is where we hold the information on the studies, series and
       %number of images in each series that we have encountered
       %a struct array with field studyUID under which there is a struct
       %array seriesData with fields seriesUID and currentImageNumber
global dataSets;

             %initialize empty data sets
dataSets = struct('studyUID','','seriesData',{});
            %initialize directory stack (push new directories onto our
            %stack and then search inside these directories) the use of a stack
            %avoids recursion
tmp = dir(inputDirectoryName);
dirIndex = find(strcmp('.',{tmp.name})==1);
dirStack = tmp(dirIndex(1));
dirStack(1).name = inputDirectoryName;
dirStackIndex = 1;

while dirStackIndex <= length(dirStack)
   currentDirectoryName = dirStack(dirStackIndex).name;
   dirStackIndex = dirStackIndex+1;
   dirData = myDir(currentDirectoryName);
   n = length(dirData);
   for i=1:n
         %an additional directory we need to check
     if dirData(i).isdir
       dirData(i).name = fullfile(currentDirectoryName, dirData(i).name);  
       dirStack(end+1) = dirData(i);
     else  %we have a file, try to read it as a DICOM file 
      try
        currentFileName = fullfile(currentDirectoryName, dirData(i).name);   
        dcmInfo = dicominfo(currentFileName);
            %add to our data structure and get the new file name. if
            %this is the first file in a series the newDirectory flag is
            %set to one.
        [newFileName, newDirectory] = addFile(dcmInfo);
        newFileName = fullfile(outputDirectoryName, newFileName);
        if newDirectory
          [newDirectoryName, fileName, fielExtension] = fileparts(newFileName);
          mkdir(newDirectoryName);
        end
        copyfile(currentFileName,newFileName);
             %show that the function is slow but alive
        fprintf('.');
      catch %we get here because we tried to get the DICOM info from a non DICOM file, just ignore it
      end
     end
   end
end

%
% Same functionality as MATLAB's 'dir' except that it omits the '.' and
% '..' directories. I don't want to assume that these are the first two
% entries even though they probably are.
%
function dirData = myDir(directoryName)

tmp = dir(directoryName);
n = length(tmp);
      %we have a non-empty directory, allocate the memory for the struct array
if n~=2      
  dirData(n-2,1) = struct('name',[],'date',[],'bytes',[],'isdir',[],'datenum',[]);
  j=1;
  for i=1:n
   if ~(strcmp(tmp(i).name,'.') || strcmp(tmp(i).name,'..'))   
     dirData(j) = tmp(i);
     j = j+1;
   end
  end
else %empty directory
  dirData = [];
end 

%
% We have a valid DICOM file, update the study/series/currentImageNumber
% data structure. If this is a new series then the new directory flag is
% set to one, otherwise it is zero.
%
function [newFileName, newDirectory] = addFile(dcmInfo)
  global dataSets;

  currentStudyUID = dcmInfo.StudyInstanceUID;
  currentSeriesUID = dcmInfo.SeriesInstanceUID;
  newDirectory = 0; 
  studyIndex = 0;
  seriesIndex = 0;
      
         %find the study and series
  if length(dataSets) > 0
     index = find(strcmp(currentStudyUID, {dataSets.studyUID})==1);
                  %we found the study, now find the series
     if size(index,2) == 1
       studyIndex = index(1);
       index = find(strcmp(currentSeriesUID,{dataSets(studyIndex).seriesData.seriesUID})==1);
       if size(index,2) == 1
         seriesIndex = index(1);
       end
     end
  end
      
       %a new study and series
  if studyIndex == 0
     dataSets(end+1) = struct('studyUID',currentStudyUID,...
                             'seriesData', struct('seriesUID',{currentSeriesUID},'currentImageNumber',0));
     studyIndex = length(dataSets);
     seriesIndex = 1;
     newDirectory = 1;
        %a new series in an existing study
  elseif seriesIndex == 0
    dataSets(studyIndex).seriesData(end+1).seriesUID = currentSeriesUID;
    dataSets(studyIndex).seriesData(end).currentImageNumber = 0;
    seriesIndex = length(dataSets(studyIndex).seriesData); 
    newDirectory = 1;
  end
  
  dataSets(studyIndex).seriesData(seriesIndex).currentImageNumber =...
    dataSets(studyIndex).seriesData(seriesIndex).currentImageNumber + 1;
  newFileName = fullfile(num2str(studyIndex),num2str(seriesIndex),...
                         strcat(num2str(dataSets(studyIndex).seriesData(seriesIndex).currentImageNumber),'.dcm'));    
 
