What is wrong in the code below?

I was trying to load DICOM images and write them as a single 3D array using the code below. But when I visualize it I see a band of differnent color somwhere in the middle of the 3D array, see below. If this was written correctly we do not suppose to see a band of dark color in this image volume. Where I made a mistake in my code? Any help would be appreciated.

close all; clear; clc

%%
folder = '/dicom/files';
files = dir(folder + "/*.dcm"); 

for k = 1:length(files)
    hdr = dicominfo(folder + "/" + files(k).name);
    
     if isfield(hdr, 'SliceLocation')
       slLoc(k,1) = hdr.SliceLocation;
     end
   
    imVol(:,:,k) = dicomread(folder + "/" + files(k).name);
end

%slLoc = sort(slLoc);
[slLoc,index] = sort(slLoc);                
slLoc = single(slLoc);        
imVol = single(imVol(:,:,index));           

volshow(imVol)
%volumeViewer(imVol)

[Walter Roberson]

I do not understand the double sort()?

[blues]

I think I forgot to erase the upper line (slLoc = sort(slLoc);) when copy-pasting the code. Sorry for the confusion. Even without it is not working. Not sure where I made mistake? I have also edited my question.

[Cris LaPierre]

I would find it easier to help if you could share your images. If possible, can you zip them and attach them to your post using the paperclip icon?

Answer 1

[Image Analyst]

Looks like some kind of rendering. If the middle slices are not as large as the upper and lower slices, then maybe they render as darker because the volume is dented in. Or maybe they just have values that are lower than the outer values of the upper and lower slices. I'm not sure how rendering works - is the voxel invisible if the value is zero?

By the way, I'd code it up like this, which will make it faster because it preallocates space for the volumetric image. And it uses the more robust fullfile.

close all; 
clear; 
clc

folder = '/dicom/files';
files = dir(fullfile(folder, "/*.dcm")); 
numberOfFiles = numel(files)
% Preallocate space to speed it up.
sliceLocation = zeros(numberOfFiles, 1); % If you want to store the slice locations for some reason.
volumetricImage = zeros(rows, columns, numberOfFiles, 'uint16'); % Use actual, known correct values for rows and columns and class.
% Read in all files inserting them into the proper location in the volume.
for k = 1 : numberOfFiles
    % Get this filename.
    thisFileName = fullfile(files(k).folder, files(k).name);
    fprintf('Reading %s.\n', thisFileName);
    thisImage = dicomread(thisFileName);
    
    % Optionally display this slice image.  (This will slow down the loop though.)
    imshow(thisImage, []);
    drawnow;

    % Check header for the slice location since filenames may not be read
    % in the actual numerical slice location order.
    hdr = dicominfo(thisFileName);    
    if isfield(hdr, 'SliceLocation')
        sliceLocation(k) = hdr.SliceLocation; % Assume it's an integer starting at 1???
    else
        sliceLocation(k) = k; % Hopefully this won't happen.  But just in case the slice location is not in the header.
    end
    % Insert this image into the proper plane/slice of the volumetric image.
    volumetricImage(:, :, sliceLocation(k)) = thisImage;
end
% Do we really need to cast to single?
volumetricImage = single(volumetricImage); 
% Render volume and display it.
volshow(imVol)
%volumeViewer(imVol)

[Image Analyst]

Sorry, no.

Edit: Answer is to the following question

Thank you Image Analysis for a better version of code. As I was working with PET iamges - I think I have missed the Rescale Slope information to incorporate it in the image volume. And the Rescale Intercept for the PET images are zero. SO, I just need to multiply individual slices by its corresponding Rescale Slope values. The Rescale Slope values seems to be slice specific. How can I do this? As I listed slice locations - I have listed Rescale Slope as: if isfield(hdr, 'RescaleSlope') ResSlope(k) = hdr.RescaleSlope; else ResSlope(k) = k; end Now my code looks like this: But I am unsure how to use Rescale Slope!

close all; 
clear; 
clc 
%% folder = '/dicom/files'; 
files = dir(fullfile(folder, "/*.dcm")); 
numberOfFiles = numel(files); 
slLoc = zeros(numberOfFiles, 1); 
imVol = zeros(256, 256, numberOfFiles, 'uint16'); 
for k = 1 : numberOfFiles 
    thisFileName = fullfile(files(k).folder, files(k).name); 
    %fprintf('Reading %s.\n', thisFileName); 
    thisImage =
    ...
end