% Example 2 - Isolate Touching Object
% -----------------------------------
% start with clean slate
close all, clear all, clc
% read image (display upper left)
afm = imread('afmsurf.tif');
figure('pos',[13 70 672 609])
subplot(331), imshow(afm), title('afm surface map')
shg, pause
% create disk-shaped structure element
se = strel('disk', 15);
% perform bottom hat filter
Ibot = imbothat(afm, se);
subplot(332), imshow(Ibot, []), title('bottom hat')
shg, pause
% perform top hat filter
Itop = imtophat(afm, se);
subplot(333), imshow(Itop, []), title('top hat')
shg, pause
% subtract images to enhance contrast (lower right)
Ienhance = imsubtract(imadd(Itop, afm), Ibot);
subplot(334), imshow(Ienhance), title('difference')
shg, pause
% convert objects of interest
Iec = imcomplement(Ienhance);
subplot(335), imshow(Iec), title('complement')
shg, pause
% detect imtensity valleys
Iemin = imextendedmin(Iec, 22);
subplot(336), imshow(Iemin), title('extended minima')
shg, pause
Iimpose = imimposemin(Iec, Iemin);
subplot(337), imshow(Iimpose), title('imposed minima')
shg, pause
% watershed segmentation
wat = watershed(Iimpose);
subplot(338), imshow(wat), title('watershed')
shg, pause
% label separate regions
rgb = label2rgb(wat);
subplot(339), imshow(rgb), title('labeled')
shg, pause
% extract features from label matrix
stats = regionprops(wat, 'Area', 'Orientation');
area = [stats(:).Area];
orient = [stats(:).Orientation];
figure('pos',[54 313 658 372])
plot(area, orient, 'b*')
title('Relationship of Particle Orientation to Area')
xlabel('Area (pixels)'), ylabel('Orientation (deg)')
shg, pause
% asses shape by exploring aspect ratio of detected features
stats = regionprops(wat,'MajorAxisLength','MinorAxisLength');
figure, plot([stats.MajorAxisLength],[stats.MinorAxisLength],'.')
xlabel('Major Axis Length (pixels)')
ylabel('Minor Axis Length (pixels)')
% Copyright 2003-2010 RBemis The MathWorks, Inc.
