quantify topology of a complex surface

Hello,

I want to quantify the topology of this complex surface. I want to know the betti number of this surface. However, the topoligica data analysis learning code does not perform on large data sets. My X1, Y1 and Zplane are 576x768 doubles. Is there some toolbox or code to quantify the topology of this surface.

Thank you.

figure(1)
A = mesh(X1,Y1,Zplane);
xlim([0 84.83])
ylim([0 63.60])
xlabel('x [\mum]')
ylabel('y [\mum]')
zlabel('z [\mum]')
title('PlaneFit')

Answer 1

[Image Analyst]

There are quite a few ways to quantify a surface, like Sa, Sq, lacunarity, StdDev, brisque, niqe, piqe, etc. But I've never heard of betti number. Do you have an equation for it?

If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Chiel van Wanrooij]

The Betti numbers of a surface can be used to extract topological information about the surface, such as the number of connected components, holes, and higher-dimensional voids. The mathematical defenition can be found on wikipedia. However, this is a quite complex defenition to put in the comment here.

I did also find this: 'Once you have a triangulated surface, you can use the betti function from the Topological Data Analysis Toolbox in Matlab to compute the Betti numbers of the surface.'

[Image Analyst]

OK, then make sure you have that toolbox and then call betti. It does not seem to be a Mathworks toolbox though.

Answer 2

[Image Analyst]

Not sure if this is of any use or is what you're after but if I code up this definition from Wikipedia using gray levels, "tells us the maximum number of cuts that can be made before separating a surface into two pieces", I get this:

% Demo by Image Analyst
% Finds betti numbers for an image which is defined by Wikiepdia at
% "tells us the maximum number of cuts that can be made before separating a surface into two pieces"
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 = 22;
markerSize = 40;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = [];
baseFileName = 'cameraman.tif';
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
	% 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
grayImage = imread(fullFileName);
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	fprintf('It is not really gray scale like we expected - it is color\n');
	% Extract the blue channel.
	grayImage = rgb2gray(grayImage);
end
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Update the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage)
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
drawnow;

%--------------------------------------------------------------------------------------------------------
% Show the histogram
subplot(2, 2, 2);
imhist(grayImage);
title('Histogram of Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
grid on;
hold on;

%--------------------------------------------------------------------------------------------------------
% Get all the gray levels in the image
gls = unique(grayImage);
bettiNumbers = [];
hl = [];

for k = 2 : length(gls) - 1
	%--------------------------------------------------------------------------------------------------------
	% Create a mask
	threshold = gls(k);
% 	% Draw red line on the histogram at the current threshold location.
% 	subplot(2, 2, 2);
% 	delete(hl); % Clear last line.
% 	hl = xline(threshold, 'Color', 'r', 'LineWidth', 2)
	mask = grayImage >= threshold;
	%--------------------------------------------------------------------------------------------------------
	% Label each blob with 8-connectivity, so we can make measurements of it
	[labeledImage, numberOfBlobs] = bwlabel(mask, 8);
	% Apply a variety of pseudo-colors to the regions.
	coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
	% Display the pseudo-colored image.
	subplot(2, 2, 3);
	cla;
	imshow(coloredLabelsImage);
	% 	imshow(mask, []);
	impixelinfo;
	axis('on', 'image');
	caption = sprintf('Number of Blobs = %d', numberOfBlobs);
	title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
	drawnow;

	% If there are exactly two blobs, save this threshold as a betti number.
	if numberOfBlobs == 2
		bettiNumbers = [bettiNumbers, threshold];
	end
end
% Show the result in the command window.
bettiNumbers
msgbox('Done!')