getmesliceorder.m

function so = getmesliceorder(directory)
% function to check the slice order in a batch of slices
% PS 2018

% Private_0027_1041 = SliceLocation
% ImagePositionPatient: delicate

%directory = '/run/media/spectro/DATAPART1/MICHA_TEST/FOK0014/e217/s1542/raw';
spion = dir(directory);
filename = [directory, '/',spion(3).name]; % should be the first i-mage after files "." and ".."
hdr = spm_dicom_headers(filename);
slicespervolume = hdr{1}.ImagesInAcquisition; % get the number of slices per volume
N=4;
for j=3 : 3 + (slicespervolume*N)-1
    tmp = spion(j).name;
    filename = [directory, '/',spion(j).name];
    hdr = spm_dicom_headers(filename);
    timing (j-2) = hdr{1}.Private_0021_105e; % TimeStamp
    location (j-2) =  hdr{1}.SliceLocation;
end

figure, plot(timing(1:slicespervolume), location(1:slicespervolume),'r.'); hold on
plot(timing(slicespervolume+1:slicespervolume*2), location(slicespervolume+1:slicespervolume*2),'g.');
plot(timing(slicespervolume*2+1:slicespervolume*3), location(slicespervolume*2+1:slicespervolume*3),'b.');
title ('Slice Y-positions of three volumes');

% extract slice order of 2nd volume
vector = timing(slicespervolume+1:slicespervolume*2);
[Y,so] = sort(vector);

% Print
str = ['The slice order here is: \n', num2str(so),'.\n'];
fprintf(str);
return