What is wrong in the code below?

1 次查看（过去 30 天）
blues
blues 2022-7-22
编辑： Cris LaPierre 2022-7-26
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)
  3 个评论
显示 1更早的评论
blues
blues 2022-7-22
编辑：blues 2022-7-22
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
Cris LaPierre 2022-7-23
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Image Analyst
Image Analyst 2022-7-23
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)
  1 个评论
Image Analyst
Image Analyst 2022-7-25
编辑：Cris LaPierre 2022-7-26
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Image Processing and Computer Vision 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by