How can I improve the image segmentation and outlining of a layered object?

I'm interested in determining the roughness at interface of two distinct layered materials from a series of images. With the code attached, I have been able to segment the brightly contrasting material from the image by k-means clustering, followed by a few closing and eroding operations on the resulting binary image. Although it works quite well at tracing the interface in some cases (as seen in the good example), but in others (like in the bad example) it is not as good.

I was wondering if anyone could give me some pointers or if there was a better way to go about this? It is also possible to do by manually thresholding but I want to automate the process.

Good example of the outline (in black) tracing the interface of the layers:

Bad example:

Code:

% Convert to gray scale if needed.
[rows, columns, numberOfColorChannels] = size(image);
if numberOfColorChannels == 3
    image = rgb2gray(image);
end
%Adjust data to span data range and adjust contrast
imstrlim = stretchlim(image);
image = imadjust(image,[imstrlim(1,1) imstrlim(2,1)]);

% Select number of clusters for kmeans clustering (number of objects in the
% image)
numberOfClusters = 4;
% Do kmeans clustering on the gray scale image.
grayLevels = double(image(:)); % Convert to column vector.
[clusterIndexes, clusterCenters] = kmeans(grayLevels, numberOfClusters,...
    'distance', 'sqEuclidean','MaxIter',10000, "Replicates",3);
labeledImage = reshape(clusterIndexes, rows, columns);
% Optional - for fun and to distinguish the classes better by giving each a unique color.
% coloredLabels = label2rgb (labeledImage, 'hsv', 'k'); % pseudo random color labels
% imshow(coloredLabels)

%Sort clusters by gray level - brightest level is selected for the
%'bulk' material of interest:
bulk_cluster = find(clusterCenters==max(clusterCenters));
BW = labeledImage == bulk_cluster;

%Select bulk material
BW = bwareafilt(BW,1);

% Close mask with disk
radius = 1;
decomposition = 0;
se = strel('disk', radius, decomposition);
BW = imclose(BW, se);

% Fill holes
BW = imfill(BW, 'holes');

% Erode mask with disk
radius = 2;
decomposition = 0;
se = strel('disk', radius, decomposition);
BW = imerode(BW, se);
BW = imfill(BW, 'holes');

%Leave outline of BW image
BW_outline = bwmorph(BW,'remove');

%Comparing Original Image with Outline
figure;
imshowpair(image,BW_outline,'diff')

Answer 1

[Image Analyst]

@NATHAN SUTEMIRE Not sure how I missed this 3 years ago, but here is my solution (better late than never). It uses kmeans to find the different brightness ranges and then takes the brightest one (as you already had). But then it does not use morphological operations, which will just smooth out the boundary and make it less accurate. Then it finds the top surface of the bulk material. Then to find the variation in it, it fits the y values to a line (to accomodate possible tilt in the sample) and then gets the distance from the actual surface to the fitted line. Then it computes a roughness index by taking the standard deviation of the distances. You can substitute a different roughness metric, such as rms or mad, if you want. I used your "rawImage_bad" image. Note how the red boundary line now follows it more accurately. How accurate it is, is dependent on kmeans. So if you trust the kmeans classification result, it will be 100% accurate.

% Demo by Image Analyst
% Initialization steps:
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%--------------------------------------------------------------------------------------------------------
%    READ IN TEST IMAGE
folder = pwd;
baseFileName = "rawImage_bad.png";
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
% Read in image file.
grayImage = imread(fullFileName);
% Get size
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels == 3
    % Convert to gray scale if needed.
    grayImage = rgb2gray(grayImage);
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
axis('on', 'image');
impixelinfo;
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
g.Name = 'Demo by Image Analyst';
g.NumberTitle = 'off';
drawnow;

%--------------------------------------------------------------------------------------------------------
% Use kmeans to find the brightest region.
% Select number of clusters for kmeans clustering 
% (number of different brightness ranges in the image).
numberOfClusters = 4;
% Do kmeans clustering on the gray scale image.
grayLevels = double(grayImage(:)); % Convert to column vector.
[clusterIndexes, clusterCenters] = kmeans(grayLevels, numberOfClusters,...
    'distance', 'sqEuclidean','MaxIter',10000, "Replicates",3);
labeledImage = reshape(clusterIndexes, rows, columns);
% Optional - for fun and to distinguish the classes better by giving each a unique color.
coloredLabels = label2rgb (labeledImage, 'hsv', 'k'); % pseudo random color labels
subplot(2, 2, 2);
imshow(coloredLabels);
axis('on', 'image');
impixelinfo;
title('Classes Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Sort clusters by gray level - brightest level is selected for the
% 'bulk' material of interest:
bulk_cluster = find(clusterCenters == max(clusterCenters));
mask = labeledImage == bulk_cluster;

%--------------------------------------------------------------------------------------------------------
% Take largest blob only.
mask = bwareafilt(mask,1);
% Fill any holes.
mask = imfill(mask, 'holes');
% Display the mask image.
subplot(2, 2, 3);
imshow(mask);
axis('on', 'image');
impixelinfo;
title('Bulk Material Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
% Use bwboundaries to get the boundaries.  This is better than just scanning the top surface because
% some "peninsulas" may curve over and have an "underneath" "bay".
% Here is where we actually get the boundaries for each blob.
boundaries = bwboundaries(mask);
% boundaries is a cell array - one cell for each blob.
% In each cell is an N-by-2 list of coordinates in a (row, column) format.  Note: NOT (x,y).
% Column 1 is rows, or y.    Column 2 is columns, or x.
numberOfBoundaries = size(boundaries, 1); % Count the boundaries so we can use it in our for loop
rc = boundaries{1}; % Extract (x,y) coordinates from cell.  It's in (row, column) order.
xy = flip(rc, 2); % Convert from (row, column), which is (y, x) to (x, y) order.
% Get rid of sides and bottoms
bottomIndexes = xy(:, 2) == rows; % Find which indexes are along the bottom of the image.
xy(bottomIndexes, :) = []; % Delete coordinates that are along the bottom of the image.
leftIndexes = xy(:, 1) == 1; % Find which indexes are along the left side of the image.
xy(leftIndexes, :) = []; % Delete coordinates that are along the left side of the image.
rightIndexes = xy(:, 1) == columns; % Find which indexes are along the right side of the image.
xy(rightIndexes, :) = []; % Delete coordinates that are along the right side of the image.

%--------------------------------------------------------------------------------------------------------
% Display the image again.
subplot(2, 2, 4);
imshow(grayImage, []);
impixelinfo;
axis('on', 'image');
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
% Plot the top edge over the image in red.
hold on;
x = xy(:, 1);
y = xy(:, 2);
plot(x, y, 'r-');
hold off;

%--------------------------------------------------------------------------------------------------------
% Find the variation of the top curve by fitting the data to a line
% (to accomodate possible tilt of the surface) and then finding the
% residuals of the vertical distance from the fitted line.  
% Then take the standard deviation of the residuals.
coefficients = polyfit(x, y, 1);
fittedY = polyval(coefficients, x);
residuals = y - fittedY;
roughnessIndex = std(residuals)
caption = sprintf('Top Edge shown in red.  Roughness Index = %f.', roughnessIndex);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');