3

One of the most widely used homogeneity measurements is the Shannon entropy: enter image description here

where p is the normalized histogram of an image of L grey levels.

We can measure such entropy by using not only image intensities but also image local gradients, since homogeneous images not only exhibit well ordered intensities but also well clustered very low gradient values in homogeneous regions ( J.V. Manjon et al, “A Nonparametric MRI Inhomogeneity Correction Method”, Medical Image Analysis, 2007).

Let Y be an image with M pixels and L1 grey levels, and let G be the associated image corresponding to the magnitude of the local gradient with L2 grey levels. The intensity gradient joint histogram is defined as: enter image description here

where δ is the Kronecker delta function. So the normalized intensity-gradient joint histogram: enter image description here

Therefore the entropy associated to the intensity-gradient joint histogram is: enter image description here

I need to calculate the above entropy for biomedical image data: https://i.stack.imgur.com/I4hf4.png. I found this discussion useful: Mutual information and joint entropy of two images - MATLAB, but I don't know if the joint histogram and the joint entropy calculated (by @rayryeng) in this discussion are the same as what I need, and if not, how I could calculate this entropy using Matlab? Any help is as always appreciated.

Community
  • 1
  • 1
margol
  • 259
  • 4
  • 12
  • Another possibly useful reference: http://stackoverflow.com/q/26512778/2778484. – chappjc Jan 02 '15 at 17:54
  • @chappjc - Wow... that code looks an awful lot like mine lol. They just changed the variable names. They even have the redundant `ones` call in `accumarray` which isn't needed. I did that in my first attempt in that link that the OP linked us to. I love how people try to post code and claim that it's their's. – rayryeng Jan 02 '15 at 17:58
  • Haha, exactly! That's how it goes when you post on SO! I was trying to be funny, but possibly my humor was a little too dry. But perhaps my comment there provides a couple alternatives to the joint entropy computation. – chappjc Jan 02 '15 at 18:00
  • @chappjc - oh I got the joke lol. I was trying to be funny back! Speaking of which, this has given me the idea to update that post to handle floating-point data. Thanks for nudge in the right direction, and happy new year! – rayryeng Jan 02 '15 at 18:03

1 Answers1

2

Yes, you can still use my post.

Looking at your question above, the Kronecker Delta function is used such that for each i and j in your joint histogram, you want to search for all values where we encounter an intensity i in the image as well as the gradient value j in the same spatial location. We count how many times these are encountered, and that goes into the ith row and jth column entry of the joint histogram.

The post you linked to is doing exactly the same thing, but the second image is just a regular intensity image. All you have to do is replace the second image with your gradient image so you can most certainly use my post.

The only thing you have to do is set im1 = Y and im2 = G. If you look at the post you linked to, I've just modified the post where the joint histogram is being calculated to make it more efficient. I unnecessarily declared a vector of all ones when I didn't have to.

The post you are referring to assume that both Y and G will be integer-valued. However, if your gradient values in G are not integer (which will most likely be the case), you can still use my post, but you'd have to assign every unique double value to a unique ID and use this ID array as input into accumarray. Therefore, what you'd have to do is use the third output of unique to help you facilitate this. This third output will give a unique ID for each unique floating point value that is encountered in G. Borrowing the code from my post, this is what you would do for each situation:

Integer-valued Gradient

indrow = double(Y(:)) + 1;
indcol = double(G(:)) + 1;
jointHistogram = accumarray([indrow indcol], 1);
jointProb = jointHistogram / length(indrow);
indNoZero = jointHistogram ~= 0;
jointProb1DNoZero = jointProb(indNoZero);
jointEntropy = -sum(jointProb1DNoZero.*log2(jointProb1DNoZero));

Floating-point Gradient

indrow = double(Y(:)) + 1;
[~,~,indcol] = unique(G(:)); %// Change
jointHistogram = accumarray([indrow indcol], 1);
jointProb = jointHistogram / length(indrow);
indNoZero = jointHistogram ~= 0;
jointProb1DNoZero = jointProb(indNoZero);
jointEntropy = -sum(jointProb1DNoZero.*log2(jointProb1DNoZero));

Note with the above where id would automatically be cast to double and it would automatically start at 1, so there's no need to offset by 1 like we did for our image intensities.

Good luck!

rayryeng
  • 102,964
  • 22
  • 184
  • 193
  • 1
    @margol - My pleasure as always. Good luck! – rayryeng Jan 02 '15 at 05:18
  • Thank you again Ray for all your efforts and help. – margol Jan 02 '15 at 16:29
  • Ray as I mentioned in my question I have, for example, an image: http://i.stack.imgur.com/I4hf4.png. This is Y in your code. What is G? How could I find this gradient image? – margol Jan 02 '15 at 17:06
  • Do you mean imgradient() function? – margol Jan 02 '15 at 17:33
  • Yes. That's it! Use that! – rayryeng Jan 02 '15 at 17:34
  • @margol - Play around with it. It will return a `double` type result, so you may need to quantize it to reduce the size of the joint histogram. That's all I can really offer because I've never done what you've talked about in your post. Good luck! – rayryeng Jan 02 '15 at 17:47
  • @margol - no such thing as a silly question when it comes to image processing! :) – rayryeng Jan 02 '15 at 17:56
  • Sure. I will play with until I reach a reliable result. I don't know why some of my comments, like the first one, were deleted after some seconds! very strange.. – margol Jan 02 '15 at 18:04
  • Ray I used G = imgradient(Y) to find the gradient image of Y. Your floating-point gradient code worked for me, since the gradient values in G are not integer. – margol Jan 02 '15 at 19:58