How to calculate the volume and sum counts

7 views (last 30 days)
mohd akmal masud
mohd akmal masud on 28 Dec 2021
Answered: Adam Danz on 28 Dec 2021
Hi all,
I have code below. Then the 3D view as below also. I also attached the workspace as matlab1 for this code (can get this file https://drive.google.com/file/d/1lXj_OhTiahKENeULn4f84VveRHUD0foo/view?usp=sharing).
% For binary images (png format, each pixel just have value 1 and 0.)
clc
clear all
dataSetDir = fullfile('C:\Users\Akmal\Desktop\I-131 256 28.02.2020\I-131 SPECT NEMA VALIDATION 01112019 256X256 26.09.2021 petang');
imageDir = fullfile(dataSetDir,'bnwaftersegmentationkahdijahkalai256');
imds = imageDatastore(imageDir);
for i = 1:90
subplot(9,10,i)
I = readimage(imds,i);
% binary
Is{i} = logical(I);
imshow(I)
end
% For 3D images spect
myFolder = ('C:\Users\Akmal\Desktop\I-131 256 28.02.2020\I-131 SPECT NEMA VALIDATION 01112019 256X256 26.09.2021 petang\dicomkhadijahkalai256');
filePattern = fullfile(myFolder, '*.dcm'); % Change to whatever pattern you need.
theFiles = dir(filePattern);
for K = 1 : length(theFiles)
baseFileName = theFiles(K).name;
fullFileName = fullfile(theFiles(K).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
% Now do whatever you want with this file name,
% such as reading it in as an image array with imread()
P(:,:,K) = dicomread(fullFileName);
P(:,:,K) = (P(:,:,K)) .* uint16( Is{K} );
% P(:,:,K) = double(uint8(K).*uint8(Is{K}));
end
% Get a list of all files in the folder with the desired file name pattern.
myFolder = ('D:\131label\CTKhadijahKalai');
filePattern = fullfile(myFolder, '*.dcm'); % Change to whatever pattern you need.
theFiles = dir(filePattern);
for L = 1 : length(theFiles)
baseFileName = theFiles(L).name;
fullFileName = fullfile(theFiles(L).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
% Now do whatever you want with this file name,
% such as reading it in as an image array with imread()
RZ(:,:,L) = dicomread(fullFileName);
end
RZ=double(RZ);
[N M L]=size(RZ);
rats=size(RZ)./size(P);
[x,y,z]=meshgrid(rats(1):rats(1):M,rats(2):rats(2):N,rats(3):rats(3):L);
RZ=interp3(RZ,x,y,z);
RZ = flip(RZ,3);
e=10000;%isovalue, air at 0 (Faceskin800, air1250)
figure
whitebg('black')
colordef none
axis off
pause
fprintf('\nFull Reconstruction. Please wait...\n');
% for e=min(min(min(RZ))):1000:max(max(max(RZ)))
% plot3(0.5,0.5,0.5);
% hold;
q = patch(isosurface(RZ,e));
axis equal
set(q,'FaceColor','w','EdgeColor','none');
alpha(q,0.2)
material shiny
hold
% q = patch(isosurface(RZ,20000));
% axis equal
% set(q,'FaceColor','g','EdgeColor','none');
% alpha(q,1)
% material shiny
camlight('headlight')
% title([num2str(e)]);
view([-22 -4])
% axis([0 M 0 N 0 L]);
% camlight('left')
% pause(0.01)
% grid
% end
p = patch(isosurface(smooth3(P)));
axis equal
set(p,'FaceColor','b','EdgeColor','none');
alpha(p,1)
hold
Then I used this code below to calculate the volume every lesion. Just right click at the lesion that I want, then export data cursor to the workspace,supposely I can get the volume. I also attached the workspace after run this below code as matlab2 (can get this file https://drive.google.com/file/d/1IfPDVuCUPCdjbU0c_1ueMMQeT-V_vy-M/view?usp=sharing).
[spect map]=dicomread('khadijahkalai256.dcm');
info = dicominfo('khadijahkalai256.dcm');
% gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3
% spect(spect<50)=nan;
e=10000;%isovalue, air at 0 (Faceskin800, air1250)
clf
whitebg('black')
colordef none
axis off
fprintf('\nFull Reconstruction. Please wait...\n');
% for e=min(min(min(RZ))):1000:max(max(max(RZ)))
% plot3(0.5,0.5,0.5);
% hold;
% % % % q = patch(isosurface(RZ,e));
% % % % axis equal
% % % % set(q,'FaceColor','w','EdgeColor','none');
% % % % alpha(q,0.2)
% % % % material shiny
% % % % hold
% % % % q = patch(isosurface(RZ,20000));
% % % % axis equal
% % % % set(q,'FaceColor','g','EdgeColor','none');
% % % % alpha(q,1)
% % % % material shiny
camlight('headlight')
% title([num2str(e)]);
view([-22 -4])
% axis([0 M 0 N 0 L]);
% camlight('left')
% pause(0.01)
% grid
% end
% spectmax=max(max(max(spect)));
% spectmin=min(min(min(spect)));
% spectrange=spectmax-spectmin;
p = patch(isosurface(P));
axis equal
set(p,'FaceColor','c','EdgeColor','none');
alpha(p,1)
[ymax xmax zmax]=size(P);
clc
disp(['-At this point, you can use all the figure tools to zoom, pan'])
disp([' and rotate your object as you wish.'])
disp([' '])
disp(['-To eliminate unnecessary object, please use data cursor tool '])
disp([' to select a point of ROI.'])
disp(['-Data tip of X, Y and Z appeared, right click and select '])
disp([' ''Export Cursor Data to Workspace''.'])
disp(['-Click OK on the pop-up window then click Enter on command '])
disp([' window. '])
ROI=1;
for a=1:ROI
a;
pause
init(a,1:3)=round([cursor_info.Position(2) cursor_info.Position(1) cursor_info.Position(3)]);
end
%init=[round(38/75*ymax) round(36/71*xmax) round(80/200*zmax)] ;%[y x]
[points dim]=size(init);
count=1;
disp([' ']);
disp(['-Wait until volume marching in Figure 1 is finished. '])
% disp(['Exit loop by clicking on the figure and then press Enter.'])
while count<=points
if 1<init(count,1) && init(count,1)<ymax && 1<init(count,2) ...
&& init(count,2)<xmax && 1<init(count,3) && init(count,3)<zmax
if P(init(count,1)+1,init(count,2),init(count,3))==1
init(end+1,1)=init(count,1)+1;
init(end,2)=init(count,2);
init(end,3)=init(count,3);
P(init(count,1)+1,init(count,2),init(count,3))=2;
end
if P(init(count,1)-1,init(count,2),init(count,3))==1
init(end+1,1)=init(count,1)-1;
init(end,2)=init(count,2);
init(end,3)=init(count,3);
P(init(count,1)-1,init(count,2),init(count,3))=2;
end
if P(init(count,1),init(count,2)+1,init(count,3))==1
init(end+1,1)=init(count,1);
init(end,2)=init(count,2)+1;
init(end,3)=init(count,3);
P(init(count,1),init(count,2)+1,init(count,3))=2;
end
if P(init(count,1),init(count,2)-1,init(count,3))==1
init(end+1,1)=init(count,1);
init(end,2)=init(count,2)-1;
init(end,3)=init(count,3);
P(init(count,1),init(count,2)-1,init(count,3))=2;
end
if P(init(count,1),init(count,2),init(count,3)+1)==1
init(end+1,1)=init(count,1);
init(end,2)=init(count,2);
init(end,3)=init(count,3)+1;
P(init(count,1),init(count,2),init(count,3)+1)=2;
end
if P(init(count,1),init(count,2),init(count,3)-1)==1
init(end+1,1)=init(count,1);
init(end,2)=init(count,2);
init(end,3)=init(count,3)-1;
P(init(count,1),init(count,2),init(count,3)-1)=2;
end
end
count=count+1;
[points dim]=size(init);
if round(count/20)==count/20
for b=1:points
c(init(b,1),init(b,2),init(b,3))=1;
end
clf
c(end+1,:,:)=0;
c(:,end+1,:)=0;
c(:,:,end+1)=0;
% plot3(0.5,0.5,0.5);
p = patch(isosurface(smooth3(smooth3(c)),0.3));
axis equal
set(p,'FaceColor','m','EdgeColor','none');
alpha(1)
view(3)
axis([min(init(:,2))-3 max(init(:,2))+3 min(init(:,1))-3 max(init(:,1))+3 min(init(:,3))-3 max(init(:,3))+3]);
camlight('right')
grid
drawnow
end
end
PixelCount=size(find(c==1))
pxlsz=info.SpatialResolution;
ROIvolume=gp*pxlsz.^2*PixelCount(1)
TotalCPS=0;
for b=1:points
b
c(init(b,1),init(b,2),init(b,3))=1;
CPSlist(b)=spect(init(b,1),init(b,2),init(b,3));
TotalCPS=TotalCPS+double(CPSlist(b));
end
TotalCPS
CPSlist
Volume=points*pxlsz^2*gp
but I Got Error
Unrecognized function or variable 'c'.
Error in PixCalc (line 182)
PixelCount=size(find(c==1))
ANYONE CAN HELP ME?

Accepted Answer

Adam Danz
Adam Danz on 28 Dec 2021
The variable c (bad variable name) is only defined when this condition is met: if round(count/20)==count/20. Since the error indicates that c is not recognized after that section of code, that condition was probably not met.
Variable that are referenced outside of a condition block should never be solely defined within the condition block. Either define the variable with a default value before the condition or add an 'else' statement to define the variable when the condition(s) is not met.
If you are expecting the condition to be met, you'll need to investigate why the condition is not satisfied. Debug mode will come in handy for that.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!