How to pass DICOM images to k mean clustering algorithm

I have given my code below,I need to do the k-means segmentation using .dcm image,when i run the code I didnt get any errors,but my matlab got stuck.Is there any mistakes in my code,that will make the matlab stuck?

img = dicomread('C:\Users\Desktop\test\hii.dcm');
  imshow(img, []);

    info = dicominfo('C:\Users\Desktop\test\hii.dcm');
    disp(info)

    noise_img = imnoise(img,'salt & pepper',0.02);
    imshow(noise_img, []);

    denoise_img = medfilt2(noise_img);

    figure;imshow(denoise_img, []);title('Denoised Image');

    [out1, out2, out3, out4] = dwt2(denoise_img,'haar');
    trans_img = [out1 out2;out3 out4];

    figure;imshow(trans_img,[]);title('Trans_Image');
    inv_trans_img = idwt2(out1,out2,out3,out4,'haar');

    figure;imshow(inv_trans_img,[]);title('inv_trans_img');

    k = 5;
    [Centroid,new_cluster]=kmeans_algorithm(inv_trans_img,k);

    for i_loop = 1:k
        cluster = zeros(size(inv_trans_img));
        pos = find(new_cluster==i_loop);
        cluster(pos) = new_cluster(pos);
       figure; imshow(cluster,[]);title('K-means');
       data=cluster;

    end

K means implementation code

function [Centroid,new_cluster]=kmeans_algorithm(input_image,k)

input_image=double(input_image);
new_image=input_image;
input_image=input_image(:);
min_val=min(input_image);
input_image=round(input_image-min_val+1);
length_input_image=length(input_image);
max_val=max(input_image)+1;
hist_gram=zeros(1,max_val);
hist_gram_count=zeros(1,max_val);
for i=1:length_input_image
  if(input_image(i)>0)
      hist_gram(input_image(i))=hist_gram(input_image(i))+1;
  end;
end
IDX=find(hist_gram);
hist_length=length(IDX);
Centroid=(1:k)*max_val/(k+1);
while(true)
  old_Centroid=Centroid;
  for i=1:hist_length
      new_val=abs(IDX(i)-Centroid);
      hist_val=find(new_val==min(new_val));
      hist_gram_count(IDX(i))=hist_val(1);
  end
  for i=1:k,
      loop_count=find(hist_gram_count==i);
      Centroid(i)=sum(loop_count.*hist_gram(loop_count))/sum(hist_gram(loop_count));
  end
  if(Centroid==old_Centroid) break;end;
end
length_input_image=size(new_image);
new_cluster=zeros(length_input_image);
for i=1:length_input_image(1),
  for j=1:length_input_image(2),
      new_val=abs(new_image(i,j)-Centroid);
      loop_count=find(new_val==min(new_val));
      new_cluster(i,j)=loop_count(1);
  end
end
Centroid=Centroid+min_val-1;

[Image Analyst]

Dicom has nothing to do with it. Use the debugger to find out where your code is "getting stuck": http://blogs.mathworks.com/videos/2012/07/03/debugging-in-matlab/, or else just put in a whole bunch of fprintf() statements, or both.

[Suba Suba]

Thank you,but still I couldn't found any errors.:( Please help me out this

[Walter Roberson]

Use the debugger or insert display statements to figure out which portion of the code is getting "stuck". Then report back to us indicating which section of code it is, include all relevant code, include the error message, and describe what you mean by "stuck".

[Suba Suba]

@Walter Roberson Sorry for the late repsponse,until the K means algorithm code,its working fine,when it come to K mean code.My matlab got freeze and pages become irresponsive. Even the debugger not helping me at all.I have updated my answer with my full code,please help me out this.

[Suba Suba]

Please some one help me to solve this problem

[Image Analyst]

What are the features? It looks like you have only one feature and that is the gray level and you are specifying that there must be 5 clusters. I have a demo for that if you're interested. (On my other computer at home, not here.)

[Zainab Mohamed]

Hi, can you help me to segment a dicom medical type image using kmean?

Answer 1

[Walter Roberson]

You have

while(true)
  old_Centroid=Centroid;
  for i=1:hist_length
      new_val=abs(IDX(i)-Centroid);
      hist_val=find(new_val==min(new_val));
      hist_gram_count(IDX(i))=hist_val(1);
  end
  for i=1:k,
      loop_count=find(hist_gram_count==i);
      Centroid(i)=sum(loop_count.*hist_gram(loop_count))/sum(hist_gram(loop_count));
  end
  if(Centroid==old_Centroid) break;end;
end

This is an infinite loop that is terminated only if all of the newly computed centroid values are bit-for-bit identical with the old centroid values. That kind of programming is vulnerable to difficulties with numeric round-off (see http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F), and is vulnerable to overshoot/undershoot problems where the computed values keep cycling around the "true" solution, first one side of the solution then the other side...

You should be testing for equality within a tolerance, and you should have a "fail-safe" maximum iteration check.

(I did not check to see if that is generally correct code for kmeans clustering: the above are reasons why that implementation could fail even if the code is generally correct.)

[Suba Suba]

Thank you for the reply.

loop = abs(Centroid-old_Centroid)
if(loop<tol) break;end

Now i changed the code like above.hope it's right.But how to calculate tolerance here?

[Walter Roberson]

As a matter of legibility I would prefer all(loop<tol) rather than just (loop<tol). You are working with vectors and although using if with a vector computation is well defined in MATLAB, people tend to not know the meaning. When I look at your code, my thought becomes that you are probably not taking into account that you are testing a vector -- but if you put in the all() then I, as a reader, know that you know what you are doing.

With regards to the fail-safe: change the

while(true)

to

for iteration = 1 : 1000

or some other number. Then you know that you will never end up in an infinite loop, that you will exit after at most that many iterations. This will not force your code to always do that many: the if/break logic you have will still allow the for loop to be exited earlier.

[Suba Suba]

Thank you for the valuable information.Please tell me how to calculate tolerance here?

[Walter Roberson]

The tolerance you use is up to you, but should probably be created in relative terms rather than absolute terms. See for example http://www.mathworks.com/help/matlab/ref/ismembertol.html . You could use ismembertol pretty directly:

if ismembertol(Centroid, old_Centroid, SomeRelativeToleranceGoesHere)

[Suba Suba]

Now I got the idea about tolerance, But still I got same problem.That matlab get freeze.Now the pages are irresponsive.

[Image Analyst]

As you step through your code, at what line of code does your program become non-responsive?

[Walter Roberson]

Instead of doing that kmeans yourself, why not call

Centroid = multithresh(inv_trans_img,k);

------

I changed your code to use

for iteration = 1 : 100

and

    if ismembertol(Centroid, old_Centroid); break;end;

and put in some fprintf() statements to monitor its progress. I also changed your script to allow the dcm file to be selected, since I do not have that particular dcm file (and if I did, it would be in a different directory.)

With the .dcm files that I had handy, your code does stop after a reasonable number of iterations (42 to 61). I did notice that the last cluster (the one with the highest index) was always coming out empty, which is something you should look at further.

To test much further I will need your hii.dcm file.

[Walter Roberson]

I tested with the four .dcm files provided with the Image Processing toolkit. Your code was able to find 5 distinct clusters for those, but not with the .dcm files that I happened to have sitting around.

One difficulty with kmeans algorithms is that it is possible for two clusters to merge, leaving a centroid associated with no points. In the Mathwork's kmeans() function, the action to be taken in such a case is controlled by the EmptyAction option.

Whether clusters merge or not can depend upon the initial positions.

kmeans clustering is not certain to find the best clustering for that number of clusters: it finds a clustering relative to the initial positions it is given. It is common for kmeans to be performed multiple times with different initial positions and then to either take the single bet fit or else to do some kind of census work about which cluster any given pixel belongs to.

[Suba Suba]

@Walter Roberson

The task given to me in my academic level is to segment the tumor cells form normal brain cells and to find it's area using k means algorithm,so I cannot change the method and I have uploaded my jpeg file with this comment. For this image, I can segment the tumor cells from above code with out any errors. But to find out the area of the tumor cells from that image, I need some metadata in DICOM image.So I thought to use DICOM images instead of this JPEG image and DICOM will give more accurate results as said in this link. How can I fix this and move on.

[Suba Suba]

@ Image Analyst

Page become irresponsive from the while loop of the above K mean implementation code

[Walter Roberson]

I need the .dcm file itself. You can zip it and attach the .zip .

You should never try to do scientific analysis on a JPEG file -- not unless the purpose is to investigate the properties of JPEG files rather than (for example) find the area of tumors. The exception to this rule is if the JPEG file has been saved lossless -- which is unlikely to be the case and is unlikely to remain the case (people tend to forget that any given JPEG needs to be lossless and then process it as if lossy is okay.)

[Suba Suba]

I have attached zip folder with some .dcm images. I couldn't get exact tumor images from the internet.It's difficult to get .dcm brain tumor images.So for now,I used the above images to test my code.

[Walter Roberson]

With the changes I noted above, http://www.mathworks.com/matlabcentral/answers/303423-how-to-pass-dicom-images-to-k-mean-clustering-algorithm#comment_392167, those dcm images converge in at most 30 iterations for the ones I tested.

[Suba Suba]

So what shall I do now to solve this problem?

[Walter Roberson]

I have attached my modified version. Invoke as

cluster_img

with no arguments. You will be prompted to select a .dcm file. You will see from the printout that the algorithm does converge within a short time.

[Suba Suba]

Thank you for the valuable time you spend with this code.But it's giving error like below

??? Undefined function or method 'ismembertol' for input arguments of type 'double'.

[Walter Roberson]

was introduced in R2015a. You failed to mention that you are using an older version.

For your purposes you can define

function tf = ismembertol(A,B)
  tf = all( abs(u-v) <= 1e-6 * max(abs([A(:);B(:)])) );

This is not exactly how Mathworks defines the function, but it is good enough for your purpose. Just watch out if you upgrade to R2015a or later afterwards.

[Suba Suba]

I'm really sorry for it,I'm using MATLAB version R2011a.I cannot update to newest version,because I'm doing a group project and my group members are doing with this R2011a version. And again it's giving error.

??? Undefined function or variable 'u'.

    Error in ==> ismembertol at 2
      tf = all( abs(u-v) <= 1e-6 * max(abs([A(:);B(:)])) );

Are those variable u,v,A and B are predefined? or do I need to change this function according to my variables name?

[Walter Roberson]

 function tf = ismembertol(A,B)
   tf = all( abs(A-B) <= 1e-6 * max(abs([A(:);B(:)])) );

This is a function, so the names A and B are the local names within the function that will be used to designate the values passed in. The only reason for you to change them would be to use better more descriptive names.

[Suba Suba]

You made my day.Finally the code is working fine. Thank you so much

[Suba Suba]

Again I have a small problem,I need to save the images that I display using imshow function into a folder as JPEG images.I referd this this link,but the think is I didn't save those image as .dcm file anywhere.So can I directly save those imshow images as .JPEG.

[Walter Roberson]

Look in your script. Do you see in your final for loop, the line

figure; imshow(cluster,[]);title('K-means');

That line is not really correct. The [] in the imshow() call means that the image data in cluster should be examined, and whatever the minimum value is should be made to show to up as black, and whatever the maximum value is should be made to show up as white. That is a problem because it is always the minimum and maximum for that one cluster that is used, independently of all of the other clusters, so this is not the same as you would get by showing the whole image but masking off parts of it (because in that case the range used would reflect all of the values.)

What you should be doing is, before the loop, do

minv = min(inv_trans_img(:));
maxv = max(inv_trans_img(:));

and then the line in the loop should be

figure; imshow(cluster,[minv, maxv]);title('K-means');

so that all of the images are scaled the same way.

However, if you do this, you will encounter a problem for any DICOM which has negative values in it (which definitely happens, especially for CT scans.) Notice how in your loop you initialize

cluster = zeros(size(inv_trans_img));

and then you fill in some positions, and then you count on the positions you did not fill in to show up as black. But if the minimum value is less than 0, then when you imshow(), those zeros are going to show up as gray (or even white). The way to solve that is to instead use

cluster = minv * ones(size(inv_trans_img));

Having made those fixes, you now want to save the resulting images as JPEG. But the images will not generally range from 0 to 1 or 0 to 255, and the images might have negative values in them. imwrite() of the cluster data requesting JPG output is not going to work well on those.

What will work instead is:

   scaled_cluster = uint16( (cluster - minv) ./ (maxv - minv) .* 65535);

and now you imwrite() of scaled_cluster and ask for 'BitDepth', 16 .

All of these changes that I showed work together to the correct final output.

[Suba Suba]

I have done these changes,but as a final output I'm getting black image.Not the tumor image.and i need to save those noised and denoised images also.

[Walter Roberson]

I think it is my bedtime now.

[Suba Suba]

oh sorry.take your own time

[Suba Suba]

@Walter Roberson, I already posted my full code here,the problem now I'm facing is when I try this code for the below image,it's extracting the tumor as well as the edge the image also there.I google it but failed at the end.Please show me a way to remove those edges.Please help me out this

P.S. I have attached both the source and the final output image,with this comment. Source image is the .zip file