Am I writing an image volume correctly from the Dicom files in a directory?

Hello everyone, I want to make sure whether the below scripts are good to write a image 3D volume using the DICOM images located in a linux directory:

%  List all files in the directory
filePattern = dir('/etc/Dicom_files/');                                             

% Remove . and .. directories
im_files = filePattern;
im_files = im_files(~ismember({im_files.name}, {'.','..'}));

% For all slice locations
slLoc = zeros(numel(im_files),1);

% Pre-allocate mage volume -manually read the image slice and use that info
imVol = zeros(256, 256, 71);

dictionaryIn = "dicom-dict.txt";
dicomdict("set",dictionaryIn);

for k = 1 : numel(im_files)
    baseFileName = im_files(k).name;
    fullFileName = fullfile(im_files(k).folder, baseFileName);
    %fprintf(1, 'Now reading %s\n', fullFileName);
    
    imVol(:, :, k) = dicomread(fullFileName);

    % Visualize 3D volume
    Volume = imshow(imVol(:,:,k)); colormap gray;
    
    % Read DICOM header
    hdr = dicominfo(fullFileName, "dictionary", dictionaryIn);
     
    % Pull slice locations from all slices as a column vector
    if isfield(hdr, 'SliceLocation')
       slLoc(k,1) = hdr.SliceLocation;
    end
end

imVol = single(imVol);
slLoc = single(sort(slLoc));

Thank you for your help in advance.

Answer 1

[Simon Chan]

Looks good except the last two lines.

You know the slice location may not be in a correct order and you sort it in the last line. However, the image set is not sorted in the same order and hence there exist inconsistency between the slice location and image set imVol.

Just simply output the index when using function sort and use it to sort the image set as shown in the last line below:

for k = 1 : numel(im_files)
    baseFileName = im_files(k).name;
    fullFileName = fullfile(im_files(k).folder, baseFileName);
    %fprintf(1, 'Now reading %s\n', fullFileName);
    
    imVol(:, :, k) = dicomread(fullFileName);
    % Visualize 3D volume
    Volume = imshow(imVol(:,:,k),[]);                          % <-- Also Change here, but it's up to you
    %colormap gray;                                            % Remove this line     
    
    % Read DICOM header
    hdr = dicominfo(fullFileName, "dictionary", dictionaryIn);
     
    % Pull slice locations from all slices as a column vector
    if isfield(hdr, 'SliceLocation')
       slLoc(k,1) = hdr.SliceLocation;
    end
end
[slLoc,index] = sort(slLoc);                        % Sort the slice location and keep the index
slLoc = single(slLoc);
imVol = single(imVol(:,:,index));                   % Rearrange the images so that it matches the slice Location

[blues]

Thank you Simon but when a image volume is visualized using the volumeviewer, central part of the image volume data seems compromised! Now sure I did it correctly - could you help me check this code again?

[Walter Roberson]

The code you posted elsewhere has an extra sort() that would throw things off.