Problem in segmenting MRI volumetric data using kmeans.

I have segmented 2D mri images in .jpg format using kmeans, fuzzy C means, region grow, active contour snake. Now i want to segment DICOM image, as these image are 16 bit with 65536 intensities and non of these method is working on 3D volumetric data (28 slices), only multi-threshold value can segment the image, how can i apply the kmeans algorithms on volumetric data, i am adding my code of threshold segmentation.

    lb =  5000; % min-threshold
    ub =  7500; % max-threshold
    mriAdjust1(mriAdjust1<=lb) = 0;
    mriAdjust1(mriAdjust1>=ub) = 0;
    bw1 = mriAdjust1>0;

[Image Analyst]

How many features are you using, and what are they? Like blob intensity and blob circularity or whatever.....

[Muhammad Shoaib]

clustering blob intensities.

[Image Analyst]

Have you segmented your image into separate blobs yet? And then you got all their intensities using regionprops()? What have you done so far? How are you defining your blobs?

[Muhammad Shoaib]

yet i didn't clustered because it gives error, i am attaching my code here which works fine for 2D images but not 3D volumetric data please make the correction in this code or any other code for clustering volumetric data.

    ima = double(mri);
    k = 3;
    imshow(ima(:,:,12));

   % check image

   copy=ima;         % make a copy
    ima=ima(:);       % vectorize ima
    mi=min(ima);      % deal with negative 
    ima=ima-mi+1;     % and zero values

   s=length(ima);

   % create image histogram

    m=max(ima)+1;
    h=zeros(1,m);
    hc=zeros(1,m);

    for i=1:s
    if(ima(i)>0) h(ima(i))=h(ima(i))+1;end;
    end
    ind=find(h);
    hl=length(ind);

    % initiate centroids

    mu=(1:k)*m/(k+1);

    % start process

    while(true)

    oldmu=mu;
    % current classification  

    for i=1:hl
      c=abs(ind(i)-mu);
      cc=find(c==min(c));
      hc(ind(i))=cc(1);
    end

    %recalculation of means  

    for i=1:k, 
      a=find(hc==i);
      mu(i)=sum(a.*h(a))/sum(h(a));
    end

    if(mu==oldmu) break;end;

    end

    % calculate mask
    s=size(copy);
    mask=zeros(s);
    for i=1:s(1),
    for j=1:s(2),
    c=abs(copy(i,j)-mu);
     a=find(c==min(c));  
    mask(i,j)=a(1);
    end
    end

   mu=mu+mi-1;   % recover real range
   %  imshow(mask==3,[]);