DicomReaderAndSaverForSTEPanizer.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reads DICOM-Files exported with
% p:\doc\MeVisLab-Networks\2011\ExtractAcinusAndExport.mlab and save them
% as JPG image sequences to further peruse with the STEPanizer.com
% First version: 20.06.2011, reading and displaying the DICOM-file
% 28.07.2011: Now reads ALL the DICOM-Files in the directory selected on
% Line 16 and converts them in one batch.

clc
clear all
close all

%% 04A - 2010a - finished - they're there
%% 04B - 2010a - finished - they're there
%% 04C - 2010a - finished - they're there
%% 04D - 2010c - finished - they're there
%% 04E - 2010c - finished - they're there

%% 10A - 2010a - finished - they're there
%% 10B - 2010a - finished - they're there
%% 10C - 2010c - finished - they're there
%% 10D - 2010c - finished - they're there
%% 10E - 2010c - finished - they're there

%% 21A - 2009f - finished - they're there
%% 21B - 2010a - finished - they're there
%% 21C
%% 21D - 2010c - finished - they're there
%% 21E - 2010c - finished - they're there

%% 36A - 2010b - finished - they're there
%% 36B - 2009b - finished - they're there
%% 36C - 2010a - finished - they're there
%% 36D - 2010a - finished - they're there
%% 36E - 2009f - segmentation not possible

%% 60A
%% 60B - 2010c - finished
%% 60C 
%% 60D - 2009f - finished
%% 60E - 2009f - finished 

Scalebar = 100; % micrometer
SliceDistance = 10;

tif = 0; % read only tiff-Files (if tif=0, read DICOM)

SamplePath = uigetdir('R:\','Select *Directory* of the Sample you want to convert the DICOM-Files from MeVisLab to JPG for STEPanizer');

for Disector = [0,1];
    for DisectorThickness = [3,5,7];
        if DisectorThickness>=SliceDistance
            disp(['DisectorThickness (' num2str(DisectorThickness) ') is equal or larger than Slicedistance (' num2str(SliceDistance) '), please redefine in MATLAB-File'])
            break
        end

        disp([ 'Counting DICOM-Files in ' SamplePath ]);         
        if tif == 1
            disp('Reading TIFF files')
            filelist = dir([SamplePath filesep '*.tif']);
        else
            disp('Reading DICOM files')
            filelist = dir([SamplePath filesep '*.dcm']);
        end
        NumberOfDICOMFiles = size(filelist,1);
        disp([ 'Found ' num2str(NumberOfDICOMFiles) ' DICOMs in ' SamplePath]);

        %% Iterate through all the files found in the directory and save them do
        %% JPG-slices
        for i=1:NumberOfDICOMFiles 
            disp([ 'Saving ' filelist(i).name ' to JPG-slices']);

            PathToDICOMFile = [SamplePath filesep]; % Carried over from first manual script, just leave it for the moment...
            DICOMFileName = filelist(i).name; % Carried over from first manual script

            %% Extract SampleName and Number of Acinus and make directory to save slices into
            SampleNameStartPointer = regexp(PathToDICOMFile, 'R108', 'once');
            SampleName = PathToDICOMFile(SampleNameStartPointer:end-1);

            VoxelSizeStartPointer = regexp(DICOMFileName, 'pixelsize', 'once');
            VoxelSize = 1000*str2num(DICOMFileName(VoxelSizeStartPointer+9:end-4)); % VoxelSize in micrometer
            ScaleBarLength = round(Scalebar/VoxelSize);

            VolumeStartPointer = regexp(DICOMFileName, 'volume', 'once');
            VolumeEndPointer = regexp(DICOMFileName, 'pixelsize', 'once');
            Volume = DICOMFileName(VolumeStartPointer+6:VolumeEndPointer-2);

            AcinusStartPointer = regexp(DICOMFileName, 'acinus', 'once');
            AcinusEndPointer = regexp(DICOMFileName, '.volume', 'once');
            AcinusNumber = DICOMFileName(AcinusStartPointer+6:AcinusEndPointer-1); % remove "acinus", so we can format the number nicely
            AcinusNumber = (sprintf('%02d',str2num(AcinusNumber))); % format string to number and pad with zero if necessary
            AcinusName = [ 'acinus' AcinusNumber ];
            if Disector == 0
                AcinusPath = [PathToDICOMFile AcinusName filesep 'voxelsize' ...
                    num2str(VoxelSize) '-every' num2str(SliceDistance) 'slice' ];
            elseif Disector == 1
                AcinusPath = [PathToDICOMFile AcinusName filesep 'voxelsize' num2str(VoxelSize) ...
                    '-every' num2str(SliceDistance) 'slice-DisectorThickness-' num2str(sprintf('%1.2f',DisectorThickness * VoxelSize)) ...
                    'um-or' num2str(DisectorThickness) 'slices' ];
            else
                warndlg('Please set Disector to either 0 or 1');
                break        
            end
            [status,message,messageid] = mkdir(AcinusPath);

            %% actually read File
            disp([ 'Reading File ' num2str(i) '/' num2str(NumberOfDICOMFiles) ]);

            if tif == 1
                DICOMFile = double(loadtiff([PathToDICOMFile DICOMFileName]));
            else
                DICOMFile = double(dicomread([PathToDICOMFile DICOMFileName])); % read in DICOM File as double
            end
            DICOMFile = DICOMFile / max(max(max(DICOMFile))).*255; % scale to a maximum value of 255, still in double

            %% Show slices of the DICOM-File
            slices = size(DICOMFile,4);
        %     subplotrows = 4;
        %      figure
        %          for ctr=1:(subplotrows^2)
        %             subplot(subplotrows,subplotrows,ctr)
        %             showslice = round(slices/(subplotrows^2)*ctr);
        %             imshow(DICOMFile(:,:,showslice),[]);
        %             title(['Slice ' num2str(showslice)])
        %          end
        %     pause(0.001)

            disp('---');

            %% Write out Slices to JPG images
            % figure
            SliceCounter = 0;
            for slice = 1:SliceDistance:slices
                disp(['writing file ' num2str(slice) '/' num2str(slices)])      
                SliceCounter = SliceCounter + 1;
                % Pad CurrentSlice to square size of longer side
                CurrentSlice = ones(max(size(DICOMFile(:,:,slice)))).*255; % Make square image with larger length of original DICOM file (white square)
                if Disector == 0
                    CurrentSlice(1:size(DICOMFile(:,:,slice),1),1:size(DICOMFile(:,:,slice),2)) = DICOMFile(:,:,slice); % Write slice of DICOM file to top left corner of white square of above line
                    CurrentSlice = uint8(CurrentSlice); % convert to uint8 before saving and displaying
                    % Make Scalebar
                    CurrentSlice(size(DICOMFile,1)-10-(round(ScaleBarLength/10)):size(DICOMFile,1)-10,10:10+ScaleBarLength) = 255; % draw Scalebar of length(ScaleBarLength) in the bottom left corner, with 10 times the length of the height.
                    WriteFileName = [ SampleName '-' AcinusName '_' num2str(SliceCounter) '.jpg' ];
                    imwrite(CurrentSlice,[AcinusPath filesep WriteFileName]);
                else
                    if slice+DisectorThickness<=size(DICOMFile,4) % only try to write slice if we actually can and if the next slice (for disector) is not out of bounds.
                        % write slice as "_a.jpg"
                        CurrentSlice(1:size(DICOMFile(:,:,slice),1),1:size(DICOMFile(:,:,slice),2)) = DICOMFile(:,:,slice); % Write slice of DICOM file to top left corner of white square of above line
                        CurrentSlice = uint8(CurrentSlice); % convert to uint8 before saving and displaying
                        % Make Scalebar
                        CurrentSlice(size(DICOMFile,1)-10-(round(ScaleBarLength/10)):size(DICOMFile,1)-10,10:10+ScaleBarLength) = 255; % draw Scalebar of length(ScaleBarLength) in the bottom left corner, with 10 times the length of the height.
                        WriteFileName = [ SampleName '-' AcinusName '_' num2str(SliceCounter) '_a.jpg' ];
                        imwrite(CurrentSlice,[AcinusPath filesep WriteFileName]);
                        % write successive slice as "_b.jpg"
                        CurrentSlice(1:size(DICOMFile(:,:,slice),1),1:size(DICOMFile(:,:,slice),2)) = DICOMFile(:,:,slice+DisectorThickness); % Write 'slice+DisectorThickness' of DICOM file to top left corner of white square of above line
                        CurrentSlice = uint8(CurrentSlice); % convert to uint8 before saving and displaying
                        % Make Scalebar
                        CurrentSlice(size(DICOMFile,1)-10-(round(ScaleBarLength/10)):size(DICOMFile,1)-10,10:10+ScaleBarLength) = 255; % draw Scalebar of length(ScaleBarLength) in the bottom left corner, with 10 times the length of the height.
                        WriteFileName = [ SampleName '-' AcinusName '_' num2str(SliceCounter) '_b.jpg' ];
                        imwrite(CurrentSlice,[AcinusPath filesep WriteFileName]);
                    end
                end
            end

            %% Give out some info
            disp(['I have written ' AcinusName ' with Volume ' num2str(Volume) ' to ' AcinusPath filesep SampleName '-' AcinusName '-x.jpg']);
            disp(['I have witten every ' num2str(SliceDistance) 'th slice!'])
            disp(['The scalebar on the image is ' num2str(Scalebar) ' micrometer long.'])
            clear DICOMFile
            disp('---')
            pause(0.001)
            close all

        end
    end
end
disp('Finished!')