Save an iterated 3D pixel sum as a vector? (Successful help -> I'll paypal the value of a beer!)

Hello,

I have written a program that extracts SPECT data from regions of interest (ROI) within an image set. Part of the program iterates a mask for the ROI over a 3D space to sum all the pixel values within that volume for the set of images.

What do I need to add to the code to successfully save the pixel sum as a vector and not scalar that's replaced each iteration? I've tried all the basics and keep getting different error messages.

Thank you for any help.

(This is part of a nested for loop iterated over a series of 29 images)

%Sum over heart

sumHeart=0;

heart = ImgBig.*uint16(heartmask);           

for h=heightLower(CaseStudy,1):heightUpper(CaseStudy,1)

      for w=widthLower(CaseStudy,1):widthUpper(CaseStudy,1)
          for d=depthLower(CaseStudy,1):depthUpper(CaseStudy,1)
              sumHeart = sumHeart + heart(h,w,d);
          end
      end
  end

[Evan]

To clarify, I want to save sumHeart as a vector.

[Geoff]

Are we talking about a decent beer, or a rubbish beer? =)

[Evan]

@Geoff: Sir! I am in graduate school, and have long since moved beyond the swill to which my undergrad budget once constrained me!

We are talking about the value of a fine microbrew from your favorite local tavern :)

Answer 1

[Image Analyst]

How about this:

counter = 1; % Before the loops start.

Then inside the loop, just make sumheart an array:

sumHeart(counter) = sumHeart(counter-1) + heart(h,w,d);
counter = counter + 1;

Is that what you want?

[Evan]

In theory, yes! I tried a variation of that earlier and got an error. I tried using yours to verify, and got same error:

Attempted to access sumHeart(2);

index out of bounds because

numel(sumHeart)=1.

Error in Main_test3 (line 119)

sumHeart(counter) =

sumHeart(counter) +

heart(h,w,d);

[Richard Brown]

make sure you put counter = counter + 1 *after* the update

[Image Analyst]

Make sure it's sumHeart(counter) = sumHeart(counter-1) +

heart(h,w,d); which is the prior value plus the new values. You might also do

if counter == 1

sumheart = heart(h,w,d);

else

sumHeart(counter) = sumHeart(counter-1) + heart(h,w,d);

end

counter = counter + 1;

See if that fixes it.

[Evan]

@Image Analyst and Richard Brown, thanks for your help....that removed the error, but unfortunately, now it's returning sumHeart as an array of all zeros :( I didn't change anything else...here's what's there now:

/code

counter = 1;

%Sum over heart

sumHeart = 0;

heart = ImgBig.*uint16(heartmask);

for h=heightLower(CaseStudy,1):heightUpper(CaseStudy,1)

for w=widthLower(CaseStudy,1):widthUpper(CaseStudy,1)

for d=depthLower(CaseStudy,1):depthUpper(CaseStudy,1)

if counter == 1

sumHeart = heart(h,w,d);

else

sumHeart(counter) = sumHeart(counter-1) + heart(h,w,d);

end

end

end

end

counter = counter + 1;

[Geoff]

You need to increment the counter inside one of your loops. Possibly immediately after the sum. You know... from the look of it, you could probably do this with no loops:

sumHeart = cumsum(heart(:));

You might need to permute the dimensions of heart before you did that, so things happened in the right order.

[Geoff]

The other way to do this without a counter (but doesn't work if you have initialise the array for efficiency) is to do:

sumHeart(end+1) = ...

I do that when I'm feeling lazy and speed is irrelevant.

[Image Analyst]

The counter incrementing needs to be IMMEDIATELY after the if statement, not after a bunch of for loops:

if counter == 1

sumHeart = heart(h,w,d);

else

sumHeart(counter) = sumHeart(counter-1) + heart(h,w,d);

end

counter = counter + 1;

[Evan]

Hey guys,

Couldn't get either of those approaches to work successfully. I kept getting sumHeart as a massive array of zeros, so I ended up just going the inefficient way of iterating the entire program.

It takes an hour to run, but whatever. Anyway, thanks for all your help. Send me your paypal email addresses (to [e-mail]) and I'll buy you each a beer for your help anyway.

Cheers.

[Geoff]

If it comes out as all zeros, then I'd be suspicious of the data in your 'heart' matrix... What is the original data? I assume it's ImgBig. Is that uint16? What is the value of max(uint16(heartmask(:)))?

I suspect you have values that are being rounded to zero in your initial array. And/or you may need to cast your values to double as you sum them.

Answer 2

[Geoff]

Just gonna chuck in my 2c worth...

h = heightLower(CaseStudy,1):heightUpper(CaseStudy,1);
w = widthLower(CaseStudy,1):widthUpper(CaseStudy,1);
d = depthLower(CaseStudy,1):depthUpper(CaseStudy,1);

H = permute( heart(h,w,d), [3 2 1] );
heartSum = cumsum(H(:));

[Evan]

See above, thanks for help!