How to adjust the code to work with 3D stack of images?

Hi everyone,

I am currently trying to modify the following code to work with the 3D stack of images needed for my image analysis:

 m= imread('ImageA.tif');
 o= imread('ImageB.tif');
 m=m(:);
 o=o(:);
 common=sum(m & o); 
 union=sum(m | o); 
 cm=sum(m); % the number of voxels in m
 co=sum(o); % the number of voxels in o
 Jaccard=common/union; % može i ovako Jaccard = sum(m(:) & o(:)) / sum(m(:) | o(:));
 Dice= 2*nnz(o&m)/(nnz(o) + nnz(m));

This current code above is for analysis of only one image. I need it to be for a 3D stack of images. I came up with the following code:

 %My inputs are a 3D stack of reference images named "ground" and the images that I would like to compare against the golden truth named "segm".

 m=zeros(numImages,1); %numImages is a number of images present in 3D stack, it is the same number in both "ground" and "segm"
 o=zeros(numImages,1);

 imSize = size( ground, 1 ) .* size( ground, 2 );
 m = reshape( ground, imSize, size( ground, 3 ) );
 o = reshape( segm, imSize, size( segm, 3 ) );

 common=sum(m & o); 
 common=transpose(common);
 union=sum(m | o); 
 union=transpose(union);

 cm=sum(m);
 cm=transpose(cm);
 co=sum(o);
 co=transpose(co);

 Jaccard=zeros(numImages,1);
 Dice=zeros(numImages,1);

 Jaccard=common/union;

 Dice=(2*common)/(cm+co); %this works the same as Dice for one image analysis above.

but when I compare the results from my code for 3D stack with the "one image at a time" approach, the results are not the same. I can not figure out where is the error.

Any help would be greatly appreciated.

Thanks!

[Rik]

You need to make sure your processing steps are element-wise and that you apply the sums over the correct dimensions. Am I correct to think that you have two binary arrays of size [sizeIM(1) sizeIM(2) numImages]?

[Mario]

I am sure that the error lies somewhere within those processing steps regarding applying correct dimensions (element-wise), but can not figure it out correctly.

I basically have two binary sets of images that I named "ground" and "segm" with identical number of images inside them which I named numImages.

Answer 1

[Rik]

When you have problems with dimension, you should use the debugger to go through your code line by line and check in the Workspace window if any variable have unexpected shapes.

sizeIM=[20,30];numImages=10;

%generate random images to test the code
ground=rand([sizeIM numImages])<0.3;
segm=rand([sizeIM numImages])<0.35;

common=ground & segm;
common=squeeze(sum(sum(common,1),2));
total=ground | segm;
total=squeeze(sum(sum(total,1),2));

%calculate Jaccard index and Sorensen-Dice coefficient
Jaccard=common./total;
%direct calculation:
%Dice=(2*common)./(squeeze(sum(sum(ground,1),2))+squeeze(sum(sum(segm,1),2)));
%derivated calculation:
Dice=2*Jaccard./(1+Jaccard);

[Rik]

Did this code solve your problem? If it did, please consider marking it as accepted answer. If not, feel free to comment here with your remaining issue.

[Mario]

Hi Rik, sorry for my absence, I was a bit bussy.

Acctually when I run your code I get a third set of results, totally different from a "single image" code and my code.

I am not sure where the problem lies.

[Mario]

My guess is that the problem should be in the definition of common and total.

Perhaps to define them separately for each image in ground and segm?

And then to define the Dice and Jaccard also separately for each of those images using a for loop?

Not sure though.

[Rik]

My code gives the same result as looping through a 3D array. Your code for a single image assumes the images are read as a 2D array, while the output for imread is almost always 3D: [row,col,color].

common and total are vectors with one element for each layer of ground/segm. This is the case because sum automatically converts a logical input to a double output.

[Mario]

Hm...so, what to do now?

[Rik]

If you confirm that your image was read to a 3D array, you know your method for a single slice doesn't work.

I don't really understand your looping method. I'm not 100% sure reshape is doing here what you expect it to do, nor do you need to use it, as sum and squeeze can do what you need as well. And at least you should replace the next two lines:

Jaccard=common/union;
Dice=(2*common)/(cm+co);
%replace with:
Jaccard=common./union;
Dice=(2*common)./(cm+co);

In short, I don't understand what your second method is doing, nor do I believe your first method is correct, so it doesn't surprise me that my method gives you a different result. I believe I understand what my method is doing, so I believe it is correct. Step through my code with the debugger and try to understand what each line is doing. It is of course entirely possible I overlooked something, so if you find something that surprises you, please comment back, and I'll have another look to check if that line is correct.

[Mario]

You are probably right. I will test the code by using debugger to go through it and get back to you.

Thanks for your help thus far.

[Mario]

Hi Rik, I debugged your code and it works fine.

However, I used a small batch of images and tested your code against the first method ("one image at a time" approach) and interestingly even though the results are different for each image, the mean value of Jaccard and Dice are the same for that batch of images using your approach and the first method.

They are indentical.

[Mario]

I have another similar problem, but I will start a new thread. Hope you can help me there to.

Thanks!