how can use kurtosis in 2D gray scale image segmentation?

I wrote this code to find skewness and kurtosis for a 2D gray scale image, I wondered if is it right like this? and how can I use the kurtosis output image to segment the input image into objects?

I = imread('cameraman.tif');
figure;imshow(I), title('input image');
I = double(I);
[n m] = size(I);
skewness = zeros(size(I));
kurtosis = zeros(size(I));
M = zeros(8,1);
V = zeros(8,1);
MO3 = zeros(8,1);
MO4 = zeros(8,1);
r = 1;
for i = 1:n
  for j = 1:m
      for ii=-1:1
          iwin=ii+i;
          if iwin<1;
              iwin=1;
          end
          if iwin>n;
              iwin=n;
          end
          for jj=-1:1
              jwin=jj+j;
              if jwin<1
                  jwin=1;
              end
              if jwin>m
                  jwin=m;
              end
              if r<= 8
              M(r,1) = I(iwin,jwin);
              else
                  r = 1;
              end
          end            
      end
      %calculate the mean value
      m1 = ((sum(M))+I(i,j))/9;
      %calculate the variance
      for r = 1:8
      V(r,1) = (M(r,1) - m1)^2;
      MO3(r,1) = (M(r,1) - m1)^3;
      MO4(r,1) = (M(r,1) - m1)^4;
      end
      x = (I(i,j)-m1)^2;
      deviation = ((sum(V)+x)/9)^0.5;
      x3 = (I(i,j)-m1)^3;
      x4 = (I(i,j)-m1)^4;
      %calculate the 3rd moment (skewness) high order statistic
      moment3 = (sum(MO3)+x3)/9;
      skewness(i,j) = (moment3/(deviation)^3);
      %calculate the 4rd moment (skewness) high order statistic
      moment4 = (sum(MO4)+x4)/9;
      kurtosis(i,j) = (moment4/(deviation)^4);
  end
end
figure, imagesc(skewness), title('skewness image');
figure, imagesc(kurtosis), title('kurtosis image');

thanks in advance

Answer 1

[Image Analyst]

I'm attaching my demo on computing image moments like skew and kurtosis. I'm having trouble conceptualizing the kind of image where its kurtosis would be the best thing to use for image segmentation. Why are you considering that instead of other more traditional measurements? Can you attach your image and tell us what you'd like to measure?

I'm not sure what your triple for loop is doing because you didn't add any comments. What is it doing? What is M?

[Rasha]

Thanks Image Analyst, I read your program its so brilliant, but I tried to make it to practice on making programs.

my main segment method is watershed but I want to modify it by using kurtosis, so I tried to use it where kurtosis measure peaked or flat of an image, I thought if I can use it as a condition after or inside the watershed method to detect edges.

I use the for loop to find a kurtosis for every pixel in image by values of 8-neighbors, is it wrong?

"M" is a vector I register the neighbors values on it to be sum after last neighbor to measuring the moment of the center pixel of 8-neighbors, where I tried to implement the next equation of kurtosis

 Mk = 1/n sum((xi - m)^k),     wher       i = 1:n,       n = 9
kurtosis = M4 / (M2)^2

I tried in cameraman image because it's available in MATLAB tool, this my main image

[Image Analyst]

Kurtosis is a measure of how "peaked" or "boxy" the histogram is. I'm not sure how that will help you with thresholding or separating blobs. Perhaps you can explain your reasoning.

Can you post your locally adaptive kurtosis image? Does it show any structure, like it's different around edges or textured areas or smooth areas?

[Rasha]

isn't the peaked in morphological image means the edges places? I build my suggestion on this idea!

my kurtosis image according to this program is

[Image Analyst]

No, an edge doesn't have a peaked histogram. If you had one gray level on one side and another gray level on the other side, then the histogram would be bimodal (two humps) and what the kurtosis of that would be I have no idea. I guess it might be different than in a uniform area which would have a single hump histogram. But even a single hump can have a wide range of kurtoses, so while your image may have some unintended ad hoc effect, I just don't think the theory justifies it in general. Anyway there are better edge detectors than that, like Canny, so why not use them?

Here are some different shapes for kurtosis: https://www.mvpprograms.com/help/mvpstats/distributions/SkewnessKurtosis You can see how kurtosis describes how peaked versus "squared off"/boxy the histogram is.

[Rasha]

thank you Image Analyst for the link,

what I understand is that the kurtosis indicate a peaked for hall image by using all values of each column to produce one value per column, but I don't understand what is the useful of this in image segmentation?

so if I use the same idea for each pixel in image by using its 8-neighbors it isn't possible to detect if this pixel is an edge or not?

and this is result of the program in more images, it seems to be segmented, isn't it?

[Image Analyst]

There are probably lots of things you could look at that look different in uniform regions of different gray levels, and around edges. If you plotted skewness, it might look similar. I just don't know of any theory that says that how box-like or cusp-like the histogram is, is a good measure of edges or mean gray level. If you just applied a pseudocolor lookup table to the original gray scale image, it might look similar. So how are these images any better for segmentation than just the original gray scale image? What do I gain from thresholding this local kurtosis image as opposed to simply thresholding the original image? In fact if you look at your medical image, it seems that the edges might have higher kurtosis and thresholding that might get edges, but in the cameraman image it appears that the fairly uniform dark coat of the man has higher kurtosis and so thresholding would get the uniform coat rather than the edges. So it doesn't seem to behave consistently. Like I said, it might by chance work for some particular type of image (like your medical image if you want edges or cameramen images if you want dark regions) but I can't see any theory in general that says kurtosis is best for any kind of segmentation, and better than usual methods that do work consistently like Sobel, Canny, local standard deviation, entropy, etc.

[Rasha]

Image Analyst, I'm so sorry for late answering, and thanks for replying me

No, I don't want to replace it the original gray scale image, I just want to find method to modify the segmentation by using it, I'm trying but I be confused about the kurtosis distribution in the link you share, lets see the normal one why it's looks like this, while when I plot the kurtosis function it give me another shape no one of the three curves (leptokurtic, mesokurtic, platikurtic)?

as an example the code below for the medical image above produce the next curve. then what is this distribution in the link refer to? I understand the plot but I couldn't applied it on the curve of distribution types.

 I = rgb2gray(imread('t1.png'));
 k = kurtosis(double(I));
 plot(k)

[Image Analyst]

Well that's the whole point - your histogram doesn't look like any histogram that is makes sense to compute the kurtosis for. Sure, every histogram will have a 4th moment, but that doesn't mean that it means anything. Even if you did have a centered box-like histogram or a cusp-like histogram (that are histograms that the kurtosis makes sense for), I'm not sure you would notice any difference in the image or that the local kurtosis would be anything worth thresholding on in theory.

[Rasha]

But the same in cameraman kurtosis plot result, the curve be like this!

and I still want to know what the distribution refer to? is it refer to relation between all pixels in an image Regardless of where those pixels locations? or the relation between each pixel and the pixel after?

[Image Analyst]

So? If you run your local kurtosis filter on an image, yes, it will have a histogram. No surprise there.

The histogram is a plot of the number of pixels having the gray level (regardless of where they're located) versus the gray level. It doesn't have anything to do with the relationship (gray level) between a pixel and an adjacent pixel. That is called the gray level cooccurrence matrix (GLCM) and is done by the function graycomatrix. GLCM is often used in texture analysis.

[Rasha]

but this last two results are for the kurtosis matlab function, I wondered why it isn't look like any of distribution shape explained in link, so you see that result is normal.

backing to my program, when I use the threshold it gives bad result, so I thinking if I can found a relation between each adjacent pixels which already in kurtosis to find if they belonging each other or not, or can find steepest descent for example. do you think that is possible?

I've the privilege of discussion with you Image Analyst, and thanks too much for all your answers.

[Image Analyst]

I don't have the kurtosis function. Is that in the Stats toolbox? It probably gives the kurtosis column-by-column. So there is absolutely no reason at all to look at the value of kurtosis as a function of column number and expect it to look like a histogram. I mean, if you took the mean gray level and plotted it vs. column number would you think it should look like the gray level histogram? No, of course not because they are two totally different things.

If you want to use edge detection I recommend you use Canny or Sobel. But I'm not really even sure what you want to measure in your CT image.

[Rasha]

 I apologies for late reply Image Analyst, yes it is in matlab toolbox, as in this link 
<http://www.mathworks.com/help/stats/kurtosis.html?searchHighlight=kurtosis>
and it's like you said a function of each column in image, so I tried to make the program above to find a kurtosis for each pixel but the resulted image refuse the function imhist to show the histogram of it. where is the wrong you think?
I need to understand the kurtosis very well to use it in watershed segmentation as reprocessing or in the processing, I don't want to use edge detection ordinary methods, I know they better for edges, I want to know what benefit of taking kurtosis for whole image? and for each pixel?
thanks

[Image Analyst]

You can use blockproc() to do it. See attached examples. Replace one of the functions in the demo (like the median filter) with kurtosis. I still think there is a lack of theory behind it in general, but whatever...it's your project.