Advanced Search

Automatic Thresholding with Expectation Maximization Algorithm in Matlab

/ Published in: MatLab
Save to your folder(s)
This code uses the Expectation Maximization algorithm for automatically thresholding an image.

Expand | Embed | Plain Text

Copy this code and paste it in your HTML 
  1. %read image
  2. img1 = rgb2gray(imread('dino.jpg'));%use your image here
  3. imshow(img1);
  4.  
  5. %parameter initialization for EM-Algo
  6. temp =[];
  7. phi1 = 0.5;
  8. phi2 = 0.5;
  9. u1  = 0;
  10. u2 = 50;
  11. sig1 = 2;
  12. sig2 = 3;
  13. tempimg1 = reshape(img1,1,numel(img1));
  14. tempimg1 = double(tempimg1);
  15. img = double (img1);
  16.  
  17.  
  18. for k = 1 : 3
  19.  
  20. for i = 1 : numel(img1)
  21.        %Expectation Step 
  22.        prob1 =  -0.9181-reallog(sig1)-(  (( img(i)- u1 )^2)/(2*(sig1^2))  );
  23.        prob2 =  -0.9181-reallog(sig2)-(  (( img(i)- u2 )^2)/(2*(sig2^2))  );
  24.  
  25.        probc1 = ( prob1*phi1 ) / ( (prob1*phi1)+ (prob2*phi2) );
  26.        probc2 = ( prob2*phi2 )/( (prob1*phi1) + (prob2*phi2) );
  27.  
  28.        temp = [temp ; probc1 probc2];
  29.  
  30. end
  31.        %Maximization Step
  32.  
  33.        phi1 = ( sum(temp(:,1)) )*(1/numel(img));
  34.        phi2 = (sum(temp(:,2)))*(1/numel(img));
  35.  
  36.  
  37.        u1 = ( sum((temp(:,1)').*tempimg1) )/( sum(temp(:,1)) );
  38.        u2 = ( sum((temp(:,2)').*tempimg1) )/( sum(temp(:,2)) );
  39.  
  40.        var1 =   ( sum( ( (tempimg1 - u1).^2  ).*( temp(:,1)' ) ) ) / ( sum(temp(:,1)) ); 
  41.        sig1 = sqrt(var1);
  42.        var2 =   ( sum( ( (tempimg1 - u2).^2  ).*( temp(:,2)' ) ) ) / ( sum(temp(:,2)) );
  43.        sig2 = sqrt(var2);
  44.  
  45.        temp = [];
  46.  
  47. end
  48.  
  49. %Thresholding
  50. thresh = ((u2 + u1)/2);
  51.  
  52. for m = 1 : numel(img1)
  53.     if img1(m) > thresh
  54.         img1(m) = 255;
  55.     else
  56.         img1(m) = 0;
  57.     end
  58. end
  59.  
  60. figure;
  61. imshow(img1);

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.