cluster analysis from edge detection

Hello,

I am looking to perform a cluster analysis (K-means) for an image. I've converted the original image to a binary image and performed an edge detection, which should help with the clustering, but am stuck on how to implement a cluster analysis from the binary image. Any suggestions are greatly appreciated! The code for the binary edge detection is below and the image is attached. The end goal here is to identify patterns of wood orientation.

clear
RGB=imread('DJI_0022.jpg'); %inputs image
I=rgb2gray(RGB); %convers to grayscale

figure
imshow(I)
BW1 = edge(I,'Canny',0.6);
imshow(BW1)

Answer 1

[Image Analyst]

Assuming you have a bunch of blobs that are the tree edges, I'd call regionprops() and ask for 'Orientation'. This will give you the average angle of every blob. You can then histogram the angles if you want.

props = regionprops(edgeImage, 'Orientation');
allAngles = [props.Orientation];
histogram(allAngles);
grid on;
xlabel('Angle', 'FontSize', 20);
ylabel('Count', 'FontSize', 20);

[Anna Marshall]

Thank you! I have calculated the orientations of the wood pieces, but actually am filtering out the trees because I only want to focus on the wood orientations for clustering. With that being said, is it possible to do a similar thing as what you mentioned above, but only look at the wood orientations based on the edge detection output (filtered out to remove trees). I'm curious as to the interpretation of the histogram. Is that grouping similar angles within the image? Thanks again!

[Image Analyst]

You can filter out the trees by doing a color segmentation. Find green and mask it away.

I do a histogram because that's what makes sense. You can choose the bin width if you want. Why use kmeans()? That makes no sense at all to me. Why cluster the angles into some predetermined number of angular groups, where the angles might change from image to image? I can't imagine any use case where that would be desireable and better than a histogram, which will accurately give a distribution of angles.

[Anna Marshall]

Thank you! I'm following most of what you suggested, but don't have a clear understanding of what the histogram in the code you suggested above is providing? My research goals are trying to cluster wood orientations in the image by orientation so I was thinking that some sort of cluster analysis would achieve this, although without ever having done one. I don't want to cluster the angles into some predetermined number of groups, but would rather like to have the analysis identify where those grouping patterns are.

[Image Analyst]

You can get the angles from the histogram of the orientations. You can look at various peaks in the histogram. Look at it this way: let's say the orientations were completely random from 0 to 360 so that the histogram was flat or almost flat. What would you say k should be? 2? 3? 10? Whatever you pick is arbitrary, and the groupings would be aribtrary and changing from one image to the next. Sure if you had trees only at 30, 60, and 90 degrees, you might say k=3, but what is the likelihood you'll have 3 groupings? Chances are there will be trees at in between angles, so then where would kmeans put the centroids? Maybe at 25, 64, 87 one time, and for another image at 29, 51, 89 degrees. Why would you want that? I still maintain you would want to look at how many trees there are at every angle or in every 10 degree range, depending on how many bins you told it to use. You could use findpeaks() on the histogram counts if you want to find local groupings. But at least with a histogram you have control over the angles get counted. With kmeans you don't have that control. All you can control is how many groups there are. Even if there are no groupings, if you tell kmeans to find 5 angle groups, it will be forced to find 5 whether or not those are valid groups.

[Anna Marshall]

That makes a lot of sense! Thank you for the thorough explanation! I understand what the histogram is showing now and agree with the approach you suggested, but it doesn't quite reach my objective. I can use the histogram to see the number of pieces that fall into each angle catagory, but in thinking through this more, I think what I am looking for is something that is slightly more visual (i.e. an overlay on the original image with an outline of the angle groupings) so that I can see how those groupings appear within the image. Does this seem like a possible task in matlab?

[Image Analyst]

Yes, you can use quiver to show little arrows over each tree, if that meets your needs. I haven't used quiver but I think you can color each arrow according to its angle if you want.

This may help you: http://blogs.mathworks.com/steve/2010/07/30/visualizing-regionprops-ellipse-measurements/

Also you could draw a short line at the blob's angle over the centroid into a binary image. Then us imoverlay() to overlay the lines on the original image.

[Anna Marshall]

Very neat! Quiver sounds like what I might be looking for. I'll look into it. Thanks again!

[Image Analyst]

Anna, see this code where I drew a line over the gray scale image where the line color corresponds to the angle:

% Initialization steps.  Brute force cleanup of everything currently existing to start with a clean slate.
% Just for the demo.  If you put this into your own function you'll want to get rid of the close and clear commands.
clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
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;

rgbImage = imread('DJI_0022.jpg'); %inputs image
hFig = figure;
subplot(2, 2, 1);
imshow(rgbImage)
axis('on', 'image');
impixelinfo;
title('Original RGB Image', 'FontSize', 20);
grayImage = rgb2gray(rgbImage); %convers to grayscale
subplot(2, 2, 2);
imshow(grayImage)
axis('on', 'image');
impixelinfo;
title('Gray Scale Image', 'FontSize', 20);
edgeImage = edge(grayImage, 'Canny', 0.6);
subplot(2, 2, 3);
imshow(edgeImage)
axis('on', 'image');
title('Canny Edge Image', 'FontSize', 20);
hFig.WindowState = 'maximized';
drawnow;
% Get coordinates and plot a line over them.
props = regionprops(edgeImage, 'PixelList', 'Orientation', 'Centroid');
lineLength = 100; % Whatever you want.
subplot(2, 2, 2); % Switch to the gray scale image.
hold on;
lineColorMap = jet(181); % One color for every 1 degrees.
for k = 1 : length(props)
% 	thisList = props(k).PixelList;
% 	fprintf('\nFor blob #%d of %d:\n    ', k, length(props))
% 	for k2 = 1 : size(thisList, 1)
% 		thisx = thisList(k2, 1);
% 		thisy = thisList(k2, 2);
% 		fprintf('(%d, %d), ', thisx, thisy);
% 	end
	% Get the centroid
	x = props(k).Centroid(1);
	y = props(k).Centroid(2);
	% Get the angle.
	angle = props(k).Orientation;
	% Get the line endpoints.
	x1 = x + lineLength * cosd(angle);
	x2 = x - lineLength * cosd(angle);
	y1 = y + lineLength * sind(angle);
	y2 = y - lineLength * sind(angle);
	thisColor = lineColorMap(round(angle + 90) + 1, :);
	fprintf('Drawing line from (%.1f to %.1f) to (%.1f, %.1f)\n', x1, y1, x2, y2);
	line([x1,x2], [y1,y2], 'Color', thisColor);
	if mod(k, 100) == 0
		drawnow;
	end
end

% Plot the histogram of angles.
subplot(2, 2, 4)
allAngles = [props.Orientation];
histObject = histogram(allAngles);
grid on;
xlabel('Angle', 'FontSize', 20);
ylabel('Count', 'FontSize', 20);
title('Histogram of angles', 'FontSize', 20);

[Anna Marshall]

This is too cool!! Wow! I'm excited to play around with this more, but is exactly what I was trying to figure out a way to do. Thank you!

[Image Analyst]

You're welcome. Thanks for accepting the answer. I actually thought it was pretty cool myself. However you might want to zoom in and make sure that the Canny lines actually line up well with the logs, because to me they look like a bunch of little dots. They might -- it's just worth zooming in and checking.

[Anna Marshall]

It's very cool! One thing I did notice when zooming in is that the lines don't actually match the orientations of the wood pieces. This might be either an issue with my canny edge detection or orientations from the histogram. Are the colored lines based on the canny detection or directly from region.prop?

[Anna Marshall]

Zoomed in figure for reference!

[Image Analyst]

Doesn't look like they align that well. Try using imoverlay() to put the Canny image over the gray scale image. If it still doesn't align well, then maybe something's wrong. You might try a COSFIRE filter, Hessian filter, or a Frangi filter which look for ridges. When you think about it, a ridgeline-finding filter should be better than an edge filter. I think they're are some in the File Exchange.

[Anna Marshall]

Thanks for the suggestion! I'll give that a try.

[Anna Marshall]

Interestingly, it does look like the Canny filter lines up well with the wood orientations in the image overlay. Looking back at the code, I'm wondering if the code to find the centroid angles before plotting those lines is what is creating a discrepancy in the colored lines?

[Image Analyst]

You might try getting the PixelList, like I already showed you in your other post, and then put those x,y values into polyfit to get the angle. Something like

% Get coordinates and fit to a line.
props = regionprops(edgeImage, 'PixelList', 'Centroid');
for k = 1 : length(props)
    thisList = props(k).PixelList;
    fprintf('\nGetting angle for blob #%d of %d.\n    ', k, length(props))
    thisx = thisList(k2, 1);
    thisy = thisList(k2, 2);
    coefficients = polyfit(thisx, thisy, 1); % Fit to a line
    angles(k) = atand(coefficients(1));
end

Answer 2

[Anna Marshall]

I'll give it a try! Thanks!

[Anna Marshall]

One more question for you! I modified the code slightly to include k2=, but am no longer able to get the orientation lines to show up.

% Get coordinates and fit to a line.
props = regionprops(edgeImage, 'PixelList', 'Centroid');
for k = 1 : length(props)
    thisList = props(k).PixelList;
    fprintf('\nGetting angle for blob #%d of %d.\n    ', k, length(props));
    k2 = 1 : size(thisList, 1); 
    thisx = thisList(k2, 1);
    thisy = thisList(k2, 2);
    coefficients = polyfit(thisx, thisy, 1); % Fit to a line
    angles(k) = atand(coefficients(1));
end

[Image Analyst]

I think it should be

% Get coordinates and fit to a line.
props = regionprops(edgeImage, 'PixelList', 'Centroid');
for k = 1 : length(props)
    thisList = props(k).PixelList;
    fprintf('\nGetting angle for blob #%d of %d.\n    ', k, length(props));
    thisx = thisList(:, 1);
    thisy = thisList(:, 2);
    coefficients = polyfit(thisx, thisy, 1); % Fit to a line
    angles(k) = atand(coefficients(1));
end

[Anna Marshall]

Getting rid of the k2 makes more sense! Thanks. The funny thing is, I keep running into an error in the props=regionprops line now that 'orientation' is removed. I'm using the following as code for this section, but the error is in props = regionprops(edgeImage, 'PixelList', 'Centroid ')

% Get coordinates and plot a line over them.
imshow(grayImage)
axis('on', 'image');
title('Line Grouping', 'FontSize', 20);
hFig.WindowState = 'maximized';
drawnow;
props = regionprops(edgeImage, 'PixelList', 'Centroid '); 
lineLength = 20; % Whatever you want.
hold on;
lineColorMap = jet(181); % One color for every 1 degrees.
for k = 1 : length(props)
    thisList = props(k).PixelList;
    fprintf('\nFor blob #%d of %d:\n    ', k, length(props))
    thisx = thisList(:, 1);
    thisy = thisList(:, 2);
    coefficients = polyfit(thisx, thisy, 1); % Fit to a line 
    angles(k) = atand(coefficients(1 )); 
    thisColor = lineColorMap(round(angle + 90) + 1, :);
    fprintf('Drawing line from (%.1f to %.1f) to (%.1f, %.1f)\n', x1, y1, x2, y2);
    line([x1,x2], [y1,y2], 'Color', thisColor);
    if mod(k, 100) == 0
		drawnow;
    end
end

[Image Analyst]

Anna, this code works but it looks like you're still missing quite a few edges with that filter.

% Initialization steps.  Brute force cleanup of everything currently existing to start with a clean slate.
% Just for the demo.  If you put this into your own function you'll want to get rid of the close and clear commands.
clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
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;

rgbImage = imread('DJI_0022.jpg'); %inputs image
hFig = figure;
subplot(2, 2, 1);
imshow(rgbImage)
axis('on', 'image');
impixelinfo;
title('Original RGB Image', 'FontSize', 20);
grayImage = rgb2gray(rgbImage); %convers to grayscale
subplot(2, 2, 2);
imshow(grayImage)
axis('on', 'image');
impixelinfo;
title('Gray Scale Image', 'FontSize', 20);
edgeImage = edge(grayImage, 'Canny', 0.6);
subplot(2, 2, 3);
imshow(edgeImage)
axis('on', 'image');
title('Canny Edge Image', 'FontSize', 20);
hFig.WindowState = 'maximized';
drawnow;
% Get coordinates and plot a line over them.
imshow(grayImage)
axis('on', 'image');
title('Line Grouping', 'FontSize', 20);
hFig.WindowState = 'maximized';
drawnow;
props = regionprops(edgeImage, 'PixelList', 'Centroid'); 
lineLength = 20; % Whatever you want.
hold on;
lineColorMap = jet(181); % One color for every 1 degrees.
angles = zeros(1, length(props));
for k = 1 : length(props)
    thisList = props(k).PixelList;
    fprintf('\nFor blob #%d of %d:\n    ', k, length(props))
    thisx = thisList(:, 1);
    thisy = thisList(:, 2);
    coefficients = polyfit(thisx, thisy, 1); % Fit to a line 
    angles(k) = atand(coefficients(1 )); 
	% Get the centroid
	x = props(k).Centroid(1);
	y = props(k).Centroid(2);
	% Get the line endpoints.
	x1 = x + lineLength * cosd(angles(k));
	x2 = x - lineLength * cosd(angles(k));
	y1 = y + lineLength * sind(angles(k));
	y2 = y - lineLength * sind(angles(k));
	thisColor = lineColorMap(round(angles(k) + 90) + 1, :);
	fprintf('Drawing line from (%.1f to %.1f) to (%.1f, %.1f)\n', x1, y1, x2, y2);    line([x1,x2], [y1,y2], 'Color', thisColor);
    if mod(k, 100) == 0
		drawnow;
    end
end

% Plot the histogram of angles.
subplot(2, 2, 4)
histObject = histogram(angles);
grid on;
xlabel('Angle', 'FontSize', 20);
ylabel('Count', 'FontSize', 20);
title('Histogram of angles', 'FontSize', 20);

[Anna Marshall]

Thank you!! That works great. I'm going to play around with different filtersm particularly the ridgeline-finding filters and see if there is one that might work a bit better than the canny filter. Thanks again for all the suggestions!!

[Anna Marshall]

@Image Analyst- I have another question to throw your way! I've been playing around with this code and one thing that I'm looking to try and do is add borders that group together "clusters" of the same colors aka same angles. For example, something like the attached sketch. I've tried the boundary function, but can't get it to work quite right. Do you have any ideas on what might work?