clear all close all [baseFileName, folder] = uigetfile('*.*', 'Specify an image file'); fullImageFileName = fullfile(folder, baseFileName); if ~exist(fullImageFileName, 'file') message = sprintf('This file does not exist:\n%s', fullImageFileName); uiwait(msgbox(message)); return; end 

 %% [rgbImage storedColorMap] = imread(fullImageFileName); [rows columns numberOfColorBands] = size(rgbImage); imshow(rgbImage) 

redBand = rgbImage(:, :, 1); greenBand = rgbImage(:, :, 2); blueBand = rgbImage(:, :, 3); fontSize = 16; [redObjectsMask,maskedRGBImage] = createMask(rgbImage); 

channel2Min = 0.5900 

 redObjectsMask= uint8(redObjectsMask); imshow(redObjectsMask, []); caption = sprintf('Mask of Only\nThe Red Objects'); title(caption, 'FontSize', fontSize); 

 %% keeping area larger than 100 pix smallestAcceptableArea = 100; % Get rid of small objects. Note: bwareaopen returns a logical. redObjectsMask = uint8(bwareaopen(redObjectsMask, smallestAcceptableArea)); figure % subplot(1, 3, 1); imshow(redObjectsMask, []); caption = sprintf('bwareaopen() removed objects\nsmaller than %d pixels', smallestAcceptableArea); title(caption, 'FontSize', fontSize); 

 %% Smooth the border using a morphological closing operation, imclose(). structuringElement = strel('disk', 4); redObjectsMask = imclose(redObjectsMask, structuringElement); % subplot(1, 3,2); figure imshow(redObjectsMask, []); fontSize = 16; title('Border smoothed', 'FontSize', fontSize); 

 %% Fill in any holes in the regions, since they are most likely red also. redObjectsMask = uint8(imfill(redObjectsMask, 'holes')); % redObjectsMask_1=imfill(redObjectsMask, 'holes'); % subplot(1, 3, 3); figure imshow(redObjectsMask, []); title('Regions Filled', 'FontSize', fontSize); 

%% % You can only multiply integers if they are of the same type. % (redObjectsMask is a logical array.) % We need to convert the type of redObjectsMask to the same data type as redBand. redObjectsMask = cast(redObjectsMask, class(redBand)); maskedImageR = redObjectsMask .* redBand; maskedImageG = redObjectsMask .* greenBand; maskedImageB = redObjectsMask .* blueBand; maskedRGBImage = cat(3, maskedImageR, maskedImageG, maskedImageB); % Show the masked off, original image. figure imshow(maskedRGBImage); fontSize = 13; caption = sprintf('Masked Original Image\nShowing Only the Red Objects'); title(caption, 'FontSize', fontSize); 

 %% bw = redObjectsMask; %imbinarize(I); bw_l = logical(bw) ; %% % Obtain the boundaries of all the coins in the binary image. [B,L,N] = bwboundaries(bw,'holes'); STATS = regionprops(L, 'Centroid', 'Area', 'Perimeter','Circularity', 'Image'); % we need 'BoundingBox' and 'Extent' Centroid = cat(1, STATS.Centroid); Perimeter = cat(1,STATS.Perimeter); Area = cat(1,STATS.Area); numSidesDistance = FindNumberOfVertices(STATS, L); CircleMetric = (Perimeter.^2)./(4*pi*Area); %circularity metric SquareMetric = NaN(N,1); TriangleMetric = NaN(N,1); 

 %% % Plot the boundary for one of the coins over the original image. hf = figure; imshow(rgbImage) hold on; dividingValues = PlotTheoreticalCircularity; for k=1:length(B) boundary = B{k}; stats= STATS(k); plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2) % [rx,ry,boxArea] = minboundrect( boundary(:,2), boundary(:,1)); %x and y are flipped in images % width = sqrt( sum( (rx(2)-rx(1)).^2 + (ry(2)-ry(1)).^2)); % height = sqrt( sum( (rx(2)-rx(3)).^2+ (ry(2)-ry(3)).^2)); % aspectRatio = width/height; % if aspectRatio > 1, % aspectRatio = height/width; %make aspect ratio less than unity % end % SquareMetric(k) = aspectRatio; %aspect ratio of box sides % TriangleMetric(k) = Area(k)/boxArea; %filled area vs box area % Requires MinBoundSuite % %for each boundary, fit to bounding box, and calculate some parameters % Use |reducepoly| to reduce the number of points defining the coin boundary. p = [boundary(:,2) boundary(:,1)]; tolerance = 0.08; % choose suitable tolerance p_reduced = reducepoly(p,tolerance); % p_reduced = unique(p_reduced,'row') plot(p_reduced(:,1),p_reduced(:,2),... 'color',[0 0 1],'linestyle','-','linewidth',2,... 'marker','o','markersize',5) thisBlob = stats.Image; allBwAreas = bwarea(thisBlob); allPerimeters = stats.Perimeter; bwCircularities = (4 * pi * allBwAreas) ./ allPerimeters.^2; %============================================================== % Determine the number of sizes according to the circularity % Get the circularity of this specific blob. thisCircularity = bwCircularities %sortedCircularities(blobNumber); % See which theoretical dividing value it's less than. % This will determine the number of sides it has. numSidesCircularity = find(thisCircularity < dividingValues, 1, 'first'); % Assign a string naming the shape according to the distance algorithm. if numSidesCircularity == 3 % Blob has 3 sides. theShapeCirc = 'triangle'; elseif numSidesCircularity == 4 % Blob has 4 sides. theShapeCirc = 'square'; elseif numSidesCircularity == 5 % Blob has 5 sides. theShapeCirc = 'pentagon'; elseif numSidesCircularity == 6 % Blob has 6 sides. theShapeCirc = 'hexagon'; else % Blob has 7 or more sides. theShapeCirc = 'nearly circular'; end %============================================================== % Determine the number of sizes according to the centroid-to-perimeter algorithm % Classify the shape by the centroid-to-perimeter algorithm which seems to be more accurate than the circularity algorithm. numSidesDist = numSidesDistance(k); % Assign a string naming the shape according to the distance algorithm. if numSidesDist == 3 % Blob has 3 sides. theShapeDistance = 'triangle'; elseif numSidesDist == 4 % Blob has 4 sides. theShapeDistance = 'square'; elseif numSidesDist == 5 % Blob has 5 sides. theShapeDistance = 'pentagon'; elseif numSidesDist == 6 % Blob has 6 sides. theShapeDistance = 'hexagon'; elseif numSidesDist == 8 % Blob has 6 sides. theShapeDistance = 'octagon'; elseif numSidesDist>9 % Blob has 7 or more sides. theShapeDistance = 'nearly circular'; else theShapeDistance = 'abnormal Shape'; end poly = polyshape(p_reduced); Peri = perimeter(poly); Area_l = polyarea(p_reduced(:,1),p_reduced(:,2)); if numSidesDist>=3 || (CircleMetric(k) < 1.1) pos=[ min(p_reduced(:,1))-10,min(p_reduced(:,2))-10, abs(max(p_reduced(:,1))- min(p_reduced(:,1)))+20,abs(max(p_reduced(:,2))-min(p_reduced(:,2)))+20]; rectangle('Position',pos,'EdgeColor','r','LineWidth',3) % Combined = [CircleMetric(k), SquareMetric(k), TriangleMetric(k)]; % Txt = sprintf('C=%0.3f S=%0.3f T=%0.3f', Combined(:)); text(pos(1)+10,pos(2)+pos(4)+20, theShapeDistance,'Color','red','FontSize',14); end % poly = polyshape(p_reduced) % Peri = perimeter(poly); % Area_l = polyarea(p_reduced(:,1),p_reduced(:,2)); % % if(Area_l>0.001*(rows*columns) && Peri<Area_l) % if (TriangleMetric(k) < 0.6) || (SquareMetric(k)>0.9 || (CircleMetric(k) < 1.1)) %Area>0.008*(rows*columns) && Peri<Area % if(length(unique(p_reduced,'row'))>=3) % pos=[ min(p_reduced(:,1))-10,min(p_reduced(:,2))-10, abs(max(p_reduced(:,1))- min(p_reduced(:,1)))+20,abs(max(p_reduced(:,2))-min(p_reduced(:,2)))+20]; % rectangle('Position',pos,'EdgeColor','r','LineWidth',3) % Combined = [CircleMetric(k), SquareMetric(k), TriangleMetric(k)]; % Txt = sprintf('C=%0.3f S=%0.3f T=%0.3f', Combined(:)); % text(pos(1)+10,pos(2)+pos(4)+20, theShapeDistance,'Color','red','FontSize',14); % % end % end % end end 

bwCircularities = 0.6961 thisCircularity = 0.6961 

bwCircularities = 0.3348 thisCircularity = 0.3348 

function [BW,maskedRGBImage] = createMask(RGB) %createMask Threshold RGB image using auto-generated code from colorThresholder app. % [BW,MASKEDRGBIMAGE] = createMask(RGB) thresholds image RGB using % auto-generated code from the colorThresholder app. The colorspace and % range for each channel of the colorspace were set within the app. The % segmentation mask is returned in BW, and a composite of the mask and % original RGB images is returned in maskedRGBImage. % Auto-generated by colorThresholder app on 10-May-2020 %------------------------------------------------------ % Convert RGB image to chosen color space I = rgb2hsv(RGB); % Define thresholds for channel 1 based on histogram settings channel1Min_red = 0.900; channel1Max_red = 0.051; channel1Min_yellow = 0.103; channel1Max_yellow = 0.157; % Define thresholds for channel 2 based on histogram settings channel2Min = 0.590 %0.328; channel2Max = 1.000; % Define thresholds for channel 3 based on histogram settings channel3Min = 0.49; channel3Max = 1.000; % Create mask based on chosen histogram thresholds sliderBW = ( ((I(:,:,1) >= channel1Min_red) | (I(:,:,1) <= channel1Max_red)) | ((I(:,:,1) >= channel1Min_yellow ) & (I(:,:,1) <= channel1Max_yellow)) ) & ... (I(:,:,2) >= channel2Min ) & (I(:,:,2) <= channel2Max) & ... (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max); BW = sliderBW; % Initialize output masked image based on input image. maskedRGBImage = RGB; % Set background pixels where BW is false to zero. maskedRGBImage(repmat(~BW,[1 1 3])) = 0; end function dividingValues = PlotTheoreticalCircularity() try dividingValues = []; % Initialize fontSize = 24; % For reference, compute the theoretical circularity of a bunch of regular polygons with different number of sides. %	fprintf('Number of Sides Theoretical Circularity\n'); % Define an array with the number of sides we want to compute the circularity for. numSides = 3 : 16; for k = 1 : length(numSides) thisSideLength = numSides(k); % Compute the theoretically perfect circularity, if the polygons were perfect instead of digitized. circularity(k) = ComputeTheoreticalCircularity(thisSideLength); end % Plot the theoretical circularities on the curve with a cross. %	plot(numSides, circularity, 'b+-', 'LineWidth', 2, 'MarkerSize', 20); %	grid on; %	hold on; xl = xlim(); % Get left and right x coordinates of the graph. % Plot theoretical lines in dark red. %	darkRed = [0.85, 0, 0]; %	for k = 1 : length(numSides) %	% Make theoretical line on the plot in a magenta color. %	line(xl, [circularity(k), circularity(k)], 'Color', darkRed, 'LineWidth', 2); % %	fprintf(' %d %f\n', thisSideLength, circularity(k)); %	if k < 7 % Only print text if it's not too crowded and close together. %	% Make text with the true value % %	message = sprintf('Theoretical value for %d sides = %.4f', thisSideLength, circularity(k)); % %	text(xl(1)+0.1, circularity(k) + 0.005, message, 'Color', darkRed); %	end %	end %	% Set up figure properties: %	% Enlarge figure to full screen. %	set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); %	% Get rid of tool bar and pulldown menus that are along top of figure. %	set(gcf, 'Toolbar', 'none', 'Menu', 'none'); %	% Give a name to the title bar. %	set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') %	drawnow; % Make it display immediately. % %	title('Theoretical Circularities', 'FontSize', fontSize, 'Interpreter', 'None'); %	xlabel('Number of Sides', 'FontSize', fontSize); %	ylabel('True Circularity', 'FontSize', fontSize); % Get the midpoint between one circularity and the one for the next higher number of sides. dividingValues = conv(circularity, [1, 1]/2, 'valid'); % Prepend two zeros so that we can use this array as a lookup table where we pass in % the number of sides as an index and it tells us the dividing value between % that number of sides and one more than that. % For example, right now dividingValues(1) gives us the dividing value between 3 and 4 % and dividingValues(3) gives us the dividing value between 5 and 6 (instead of between 3 and 4). dividingValues = [0, 0, dividingValues]; % Now dividingValues(3) will give us the dividing value between 3 and 4. %	% Put up red horizontal lines at the dividing values %	hold on; %	xl = xlim(); %	darkGreen = [0, 0.5, 0]; %	for k = 1 : length(numSides)-1 %	thisSideLength = numSides(k); %	thisDividingValue = dividingValues(thisSideLength); %	h = line(xl, [thisDividingValue, thisDividingValue], 'Color', darkGreen, 'LineWidth', 2, 'LineStyle', '--'); %	%	h.LineStyle = '--'; %	% For the first 6, print the dividing value at the left just above the line. %	% After 6 it would get too crowded %	if k <= 6 %	theLabel = sprintf('Dividing value = %.4f', thisDividingValue); %	text(xl(1)+0.1, thisDividingValue + 0.005, theLabel, 'Color', darkGreen); %	end %	end catch ME errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ... ME.stack(1).name, ME.stack(1).line, ME.message); fprintf(1, '%s\n', errorMessage); uiwait(warndlg(errorMessage)); end end % https://en.wikipedia.org/wiki/Regular_polygon % Which says A = (1/4) * n * s^2 * cot(pi/n) function circularity = ComputeTheoreticalCircularity(numSides) try sideLength = 1; perimeter = numSides * sideLength; area = (1/4) * numSides * sideLength^2 / tan(pi / numSides); circularity = (4 * pi * area) / perimeter ^2; catch ME errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ... ME.stack(1).name, ME.stack(1).line, ME.message); fprintf(1, '%s\n', errorMessage); uiwait(warndlg(errorMessage)); end end % Now compute the number of vertices by looking at the number of peaks in a plot of distance from centroid. function numVertices = FindNumberOfVertices(blobMeasurements, labeledImage) try numVertices = 0; % Initialize. % Get the number of blobs in the image. numRegions = length(blobMeasurements); hFig = figure; promptUser = true; % Let user see the curves. % For each blob, get its boundaries and find the distance from the centroid to each boundary point. for k = 1 : numRegions % Extract just this blob alone. thisBlob = ismember(labeledImage, k) > 0; %	if promptUser % Let user see the image. %	cla; %	imshow(thisBlob); %	end % Find the boundaries thisBoundary = bwboundaries(thisBlob); thisBoundary = cell2mat(thisBoundary); % Convert from cell to double. % Get x and y x = thisBoundary(:, 2); y = thisBoundary(:, 1); % Get the centroid xCenter = blobMeasurements(k).Centroid(1); yCenter = blobMeasurements(k).Centroid(2); % Compute distances distances = sqrt((x - xCenter).^2 + (y - yCenter).^2); %	if promptUser % Let user see the curves. %	% Plot the distances. %	plot(distances, 'b-', 'LineWidth', 3); %	grid on; %	message = sprintf('Centroid to perimeter distances for shape #%d', k); %	title(message, 'FontSize', 15); %	% Scale y axis %	yl = ylim(); %	ylim([0, yl(2)]); % Set lower limit to 0. %	end % Find the range of the peaks peakRange = max(distances) - min(distances); minPeakHeight = 0.5 * peakRange; % Find the peaks [peakValues, peakIndexes] = findpeaks(distances, 'MinPeakProminence', minPeakHeight); % Find the valueys. [valleyValues, valleyIndexes] = findpeaks(-distances, 'MinPeakProminence', minPeakHeight); numVertices(k) = max([length(peakValues), length(valleyValues)]); % Circles seem to have a ton of peaks due to the very small range and quanitization of the image. % If the number of peaks is more than 10, make it zero to indicate a circle. if numVertices(k) > 10 numVertices(k) = 0; end %	if promptUser % Let user see the curves. %	% Plot the peaks. %	hold on; %	plot(peakIndexes, distances(peakIndexes), 'r^', 'MarkerSize', 10, 'LineWidth', 2); % %	% Plot the valleys. %	hold on; %	plot(valleyIndexes, distances(valleyIndexes), 'rv', 'MarkerSize', 10, 'LineWidth', 2); % %	message = sprintf('Centroid to perimeter distances for shape #%d. Found %d peaks.', k, numVertices(k)); %	title(message, 'FontSize', 20); % %	% The figure un-maximizes each time when we call cla, so let's maximize it again. %	% Set up figure properties: %	% Enlarge figure to full screen. %	set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); %	% Get rid of tool bar and pulldown menus that are along top of figure. %	set(gcf, 'Toolbar', 'none', 'Menu', 'none'); %	% Give a name to the title bar. %	set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') % %	% Let user see this shape's distances plotted before continuing. %	promptMessage = sprintf('Do you want to Continue processing,\nor Cancel processing?'); %	titleBarCaption = 'Continue?'; %	button = questdlg(promptMessage, titleBarCaption, 'Continue', 'Cancel', 'Continue'); %	if strcmpi(button, 'Cancel') %	promptUser = false; %	end %	end end close(hFig); catch ME errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ... ME.stack(1).name, ME.stack(1).line, ME.message); fprintf(1, '%s\n', errorMessage); uiwait(warndlg(errorMessage)); end end