Moving inside 3d image stack - Volume fraction calculator

Hi all,

So I am writing a program to determine the phase volume fraction of '1' values in a 3D binary image stack. Starting from Column 1, Row 1, the program horizontally reads each value and adds it to the total. at the end of each column it determines the total volume fraction of '1's and records it in an array. At the end of each row it goes to the next row and starts again from the first column. I would ideally want it to go to the next image/slice after this step.

The issue I am having is that when I finish scanning an image, I am trying to move to the next image in the same manner i move between columns and rows, however the code does not seem to be able to do that and I seem to be stuck at the same image. Please find attached the full code with comments. Any help would greatly be appreciated.

function [x, y] = linestack(V)

%sets position counters (pc = position column, pr = position row etcettera), y = volume fraction
y=0;
pc=1;
pr=1;
ps = 1;

%sets total number of rows, columns and slices/images
[r, c, s] = size(V);

%creates grid to store results
g = zeros(r,3);
fraction = 0;

%iterations to move through images
for ps=1:s
for pr = 1:r
for pc=1:c

%calculate volume fractions if pixel =1
    if V(pr,pc) == 0;
    pc=pc+1;

      elseif V(pr,pc) == 1;
      y=y+1;
      pc=pc+1;

      end
      %calculate volume fraction and write results in grid
  fraction = y/c;

g(pr,1) = pr;
g(pr,2) = ps;
g(pr,3) = fraction;

end
%reset volume fraction variables and move to next row and first column
fraction = 0;
y=0;
pr = pr + 1;
pc = 1;    

end
%move to next image (DOES NOT WORK HERE) and reset columns and rows to first position
ps = ps+1;

pr = 1;
pc = 1;

   end
  %display grid
  disp(g) 

end

Answer 1

[Walter Roberson]

Your code can be pretty much replaced with

   fraction = squeeze(mean(V, 2));

But in terms of the code you have: each place you have

V(pr,pc)

should be

V(pr,pc,ps)

Answer 2

[Image Analyst]

If you want the total volume fraction why can't you just sum the 1's:

volumeFraction = sum(V(:))/numel(V);

If you want the volume fraction of V as a function of location, so that you have a volume fraction value for every voxel location in the image, then simply use convn() to get the sum in a moving block (cube).

Answer 3

[Sohrab Daemi]

I have tried changing the values in the if iteration to

if true
  % code
  if V(pr,pc,ps) == 0;
  pc=pc+1;

    elseif V(pr,pc,ps) == 1;
    y=y+1;
    pc=pc+1;

    end
  end

but no difference in output.

Answer 4

[Sohrab Daemi]

sorted with a while command. Thanks for the shortcuts will consider them for my next tasks!