demo1_segment_rice.m

% Example1 - Size distribution of detected objects
% ------------------------------------------------

% start with clean slate
close all, clear all, clc

% use example image
im = imread('rice.tif');
figure('pos',[6 69 330 314])  %fig1, lower 1/2, left 1/3
imshow(im,'n'), title('original image')

% look at intensity distribution
figure('pos',[246 69 330 314])  %fig2, lower 1/2, left-mid 1/3
lev = graythresh(im); 
imhist(im), title(num2str(lev*255))
ax=axis; line(lev*255*[1 1],ax(3:4),'color','r','linewidth',2)

% apply threshold to segment light grains from dark background
BW = im2bw(im,lev);
figure('pos',[490 69 330 314])  %fig3, lower 1/2, right-mid 1/3
imshow(BW,'n'), title('threshold segmented')

% label individual grains with different colors
L = bwlabel(BW);
num_grains = max(L(:));
mid_gray = [0.5 0.5 0.5];
map = [mid_gray; jet(num_grains)];
figure('pos',[690 69 330 314])  %fig4 lower 1/2, right 1/3
imshow(L,map,'n'), title('labeled regions')

% extract area values from each detected object
stats = regionprops(L)
A = [stats.Area];
figure(2), hist(A), title('size distribution (area)')
% what's wrong with distribution?

% change colormap (more intuitive)
figure(1), colormap(jet(256))

% look at edges of detected grains
P = bwperim(BW);
im2=im; im2(P)=0;
map = [0 0 0; jet(255)];
figure(1), imshow(im2,map,'n'), title('detected objects')

% characterize nonuniform background using morphology
stats = regionprops(L,'MajorAxisLength');
avg_len = median([stats.MajorAxisLength])
BG = imopen(im,strel('disk',round(avg_len)));
figure(2), imshow(BG,jet(256),'n'), title('background')

% remove background by image subtraction
im2 = imsubtract(im,BG);
figure(3), imshow(im2,jet(256),'n'), title('difference image')

% segment after background removal
BW2 = im2bw(im2,graythresh(im2));
figure(4), imshow(BW2,'n'), title('segment (w/o background)')

% extract new areas and update distribution
L2=bwlabel(BW2); stats2=regionprops(L2,'area')
A2=[stats2.Area]; 
figure(2), hist(A2), title('new size distribution')

% compare distributions
bins = 25:50:800;
n=hist(A(:),bins); n2=hist(A2(:),bins);
figure(2), plot(bins,[n; n2]), title('size distributions')
legend('original','-background')

% now remove partial grains cut off by edges
BW3 = imclearborder(BW2);
figure(3), imshow(BW3,'n'), title('border objects removed')

% redo size distribution again after border clearing
L3=bwlabel(BW3); stats3=regionprops(L3,'area');
A3=[stats3.Area]; n3=hist(A3(:),bins);
figure(2), plot(bins,[n; n2; n3]), title('size distributions')
legend('original','-background','-borders')

% also compare first order statistics
mu_sig = [mean(A) std(A); mean(A2) std(A2); mean(A3) std(A3)]

% finally, indicate detected grains on original image
P3 = bwperim(BW3);
im3 = im;
im3(P3) = 0;
map = [0 0 0; jet(255)];
figure(3), imshow(im3,map,'n'), title('final detected grains')

% Copyright 2003-2010 RBemis The MathWorks, Inc.