One of the most widely used homogeneity measurements is the Shannon entropy:

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:

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

Therefore the entropy associated to the intensity-gradient joint histogram is:

I need to calculate the above entropy for biomedical image data: https://i.sstatic.net/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.

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 i^th row and j^th 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!