how to write loop coding for dicom image.

Dear all,

i have image dicom (as attached) for SPECT machine. But in this image, there are 72 slice (frame) can see in dicominfo.

But when i want to apply adaptthresh to every image, i have to write as code below for every slice. So then i have to change the number slice every time.

%let say I want to know the slice number 38
I = dicomread('spect128x128', 'frame', 38);
%determine threshold
T = adaptthresh(I, 0.4);
%change grey to binary
BW = imbinarize(I,T);
%open image
figure
imshowpair(I, BW, 'montage')

CC = bwconncomp(BW)
[r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
T = regionprops('table', BW,'Area','Centroid')

My question is, how i want to write for the loop code for the first slice(frame) till last slice. So that i no need to change the number of slice every time, and adaptthresh can apply for all the slice.

Please help me.

Answer 1

[Image Analyst]

Try this:

for sliceIndex = 1 : numSlices
	% Read this slice.
	fprintf('Reading slice #%d...\n', sliceIndex);
	thisImage = dicomread('spect128x128', 'frame', sliceIndex);
	% Determine threshold
	threshold = adaptthresh(thisImage, 0.4);
	% Change grey to binary
	BW = imbinarize(thisImage, threshold);
	% Display images.
	imshowpair(thisImage, BW, 'montage')
	drawnow;
	
% 	CC = bwconncomp(BW);
% 	[r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
	propTables{sliceIndex} = regionprops('table', BW, 'Area', 'Centroid')
end

[mohd akmal masud]

thank you. its function well.

1) but if i want to apply this code (as below), its look like it give the separate information. how to combine the information so that it give just one info for all slice that i need.

T = regionprops('table', BW,'Area','Centroid')

2) how to develop the image figure so that i can scroll slice by slice?

Please help me image analyst.

[Image Analyst]

Since BW was created from just a single slice, T will be unique to just that slice. If you don't use an index, T will be overwritten on every iteration, so you may want to use it somehow before it vanished at the end of the loop.

Look on the Apps tab - I think there is a Volume Viewer App that you can use to scroll through slices.

[mohd akmal masud]

alright. that matter i will be eliminate.

so for the first question, can you help me? i want to the infirmation gather for all slice in one table.

i have try this code, but pop up error

T = regionprops(thisImage, 'table', BW,'Area','Centroid')

can help me?

[Image Analyst]

You forgot to tell me the error. But try this. It works. Adapt as needed.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

numSlices = 2;
for sliceIndex = 1 : numSlices
	% Read this slice.
	fprintf('Reading slice #%d...\n', sliceIndex);
% 	thisImage = dicomread('spect128x128', 'frame', sliceIndex);
	thisImage = imread('cameraman.tif');
	subplot(2, 2, 1);
	imshow(thisImage);
	impixelinfo;
	caption = sprintf('Slice #%d', sliceIndex);
	title(caption, 'FontSize', fontSize);

	% Determine threshold.  It's an image where every pixel is the threshold.
	threshold = adaptthresh(thisImage, 0.4);
	% Convert to gray levels in the range 0-255
	threshold = threshold * 255;
	subplot(2, 2, 2);
	imshow(threshold, []);
	impixelinfo;
	title('Locally Adaptive Threshold', 'FontSize', fontSize);

	% Change grey to binary
% 	BW = imbinarize(thisImage, threshold);
	BW = thisImage > threshold;
	% Display images.  Need to multiple BW by 255 to get it into the same range as thisImage.
	subplot(2, 2, 3);
% 	imshowpair(thisImage, 255 * BW, 'montage')
	imshow(BW);
	impixelinfo;
	title('Binary Image', 'FontSize', fontSize);
	drawnow;
	
	propTables{sliceIndex} = regionprops('table', BW, 'Area', 'Centroid')
end

[mohd akmal masud]

ok, let me explain.

Question No. 1.

when i run this code as below,

for sliceIndex = 38 : 55
	% Read this slice.
	fprintf('Reading slice #%d...\n', sliceIndex);
	thisImage = dicomread('I-131101', 'frame', sliceIndex);
	% Determine threshold
	threshold = adaptthresh(thisImage, 0.4);
	% Change grey to binary
	BW = imbinarize(thisImage, threshold);
	% Display images.
    figure
	imshowpair(thisImage, BW, 'montage')
	drawnow;
	
%  	CC = bwconncomp(BW)
%  	[r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
% 	propTables{sliceIndex} = regionprops('table', BW, 'Area', 'Centroid')
    
%     [r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
 T = regionprops('table', BW,'Area','Centroid')
end

it wil give the onformation separately like this

Reading slice #43...
T =
  12×2 table
    Area                  Centroid              
    ____    ____________________________________
    301     11.5382059800664     64.468438538206
    703     70.6685633001423    23.2290184921764
    262     37.2977099236641    110.595419847328
      1                   36                 106
    298     47.2281879194631    63.6778523489933
     20                54.75                69.9
     39     64.4615384615385    120.641025641026
     51     64.9607843137255    78.0392156862745
     76     75.6973684210526    58.2236842105263
    464      103.29525862069    95.3362068965517
     44     76.0454545454545    70.9318181818182
    136     89.3602941176471    74.9191176470588
Reading slice #44...
T =
  7×2 table
    Area                  Centroid              
    ____    ____________________________________
    1772     63.627539503386    64.1111738148984
     296    47.1317567567568    63.8952702702703
      24              54.625    70.2083333333333
      51    64.9411764705882    77.7647058823529
      78    75.7948717948718    58.1282051282051
      44    76.0454545454545    70.9318181818182
     142    89.3239436619718    74.9507042253521
Reading slice #45...
T =
  7×2 table
    Area                  Centroid              
    ____    ____________________________________
    1800    63.6344444444444    64.4238888888889
     306    47.1143790849673    64.6013071895425
      23    54.5652173913044     70.304347826087
      46    64.9347826086957    77.3260869565217
      79     75.746835443038    58.0886075949367
      44    76.0454545454545    70.9318181818182
     143    89.1888111888112    74.8951048951049
Reading slice #46...
T =
  8×2 table
    Area                  Centroid              
    ____    ____________________________________
    1245    62.9269076305221     85.395983935743
     545    69.5302752293578    17.0293577981651
     312    47.3910256410256    64.8141025641026
      19    54.4736842105263    70.8421052631579
      44    64.8863636363636    77.2954545454545
      79     75.746835443038    58.0886075949367
      41                  76                  71
     148    88.5878378378378    75.5810810810811

so, i want to combine all the info area and centroid in one table only. so that all the info slice is summation for all area. For example, i take the last area with centroid (88, 75), the area for slice 43 is 136, slice 44 is 142, slice 45 is 143 and slice 46 is 148, so what i want is, just in one table, the area (88,75) is summation area for all the slice.

so i no need to see one by one tble for all the slice.

Question No. 2.

when i run the code same above, the figure appear seperately, like picture below

so can you help me to combine all the picture in one figure. so that i can scroll slice by slice to the gray and BW picture. my expected figure as below.

Hope you can help me. Thank you so much.

[Image Analyst]

To combine the individual tables into one table, use vertcat(). here's a little snippet to demo it:

for k = 1 : 3
	% Get the table for this image.  (REPLACE WITH YOUR ACTUAL CODE)
	thisTable = table(1000*rand(10, 1), 1000*rand(10, 2), 'VariableNames', {'Area', 'Centroid'})
	% Append the table for this image onto the table for ALL images.
	if k == 1
		% No table yet, so make the table for the first one the master one.
		overallTable = thisTable
	else
        % Overall table exists already so append thisTable onto that with vertcat().
		overallTable = vertcat(overallTable, thisTable)
	end
end

If you want all plots on one figure, use imtile() or subplot() and don't call figure() at all.

[mohd akmal masud]

is it like this?

for sliceIndex = 43 : 46
	% Read this slice.
	fprintf('Reading slice #%d...\n', sliceIndex);
	thisImage = dicomread('I-131101', 'frame', sliceIndex);
	% Determine threshold
	threshold = adaptthresh(thisImage, 0.4);
	% Change grey to binary
	BW = imbinarize(thisImage, threshold);
	% Display images.
    figure
	imshowpair(thisImage, BW, 'montage')
	drawnow;
	
%  	CC = bwconncomp(BW)
%  	[r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
% 	propTables{sliceIndex} = regionprops('table', BW, 'Area', 'Centroid')
    
%     [r, c] = cellfun(@(x) ind2sub(size(BW), x), CC.PixelIdxList, 'UniformOutput', 0)
 T = regionprops('table', BW,'Area','Centroid')
end


for k = 1 : 3
	% Get the table for this image.  (REPLACE WITH YOUR ACTUAL CODE)
	thisTable = table(1000*rand(10, 1), 1000*rand(10, 2), 'VariableNames', {'Area', 'Centroid'})
	% Append the table for this image onto the table for ALL images.
	if k == 1
		% No table yet, so make the table for the first one the master one.
		overallTable = thisTable
	else
        % Overall table exists already so append thisTable onto that with vertcat().
		overallTable = vertcat(overallTable, thisTable)
	end
end

[Image Analyst]

No. Just one loop.

% Processes slices 43 through 46.
for sliceIndex = 43 : 46
	% Read this slice.
	fprintf('Reading slice #%d...\n', sliceIndex);
	thisImage = dicomread('I-131101', 'frame', sliceIndex);
	% Determine threshold
	fprintf('   Getting threshold image...\n');
	threshold = adaptthresh(thisImage, 0.4);
	% Change grey to binary
	fprintf('   Thresholding slice #%d\n', sliceIndex);
	BW = imbinarize(thisImage, threshold);
	% Display images.
	figure
	imshowpair(thisImage, BW, 'montage')
	drawnow;

	% Get the measurements for this one slice.
	fprintf('   Making measurements via regionprops()...\n');
	thisTable = regionprops('table', BW,'Area','Centroid');
	% Append the table for this image onto the table for ALL images.
	fprintf('   Appending measurements to table...\n');
	if k == 1
		% No table yet, so make the table for the first one the master one.
		overallTable = thisTable;
	else
		% Overall table exists already so append thisTable onto that with vertcat().
		overallTable = vertcat(overallTable, thisTable);
	end
end