MATLAB Answers

How to measure the moon crescent

101 ビュー (過去 30 日間)
Ahmad Najmuddin
Ahmad Najmuddin 2020 年 12 月 21 日
Hi, anyone can help me how to measure the moon crescent using regionprops and how to remove the unwanted object in the image?

  0 件のコメント

サインインしてコメントする。

採用された回答

Matt J
Matt J 2020 年 12 月 22 日
編集済み: Matt J 2021 年 2 月 19 日 4:43
Right now i want to measure the area of the crescent and then viscircles on it.
Using this FEX package,
load(websave('Image','https://www.mathworks.com/matlabcentral/answers/uploaded_files/467970/Image.mat'))
BW0=BW;
imshow(BW0);
BW = bwpropfilt(BW,'EquivDiameter',2);
crescentArea=bwarea(BW)
crescentArea = 1.2870e+04
[y,x]=find(bwskel(BW));
obj=circularFit([x,y].'); %EDIT
imshow(BW);
hold on
viscircles(obj.center,obj.radius);
hold off

  5 件のコメント

表示 2 件の古いコメント
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 13 日
Matt, I had try imbinarize but an error occur. Please help me. I add my coding for preprocessing at the beginning and then I add your coding to measure the crescent but error when I want to convert the black area of the crescent into white colour (imbinarize like you said). This is the coding and the original image
%upload rgb image
rgb = imread('Rabiulawal 1439.png');
[J,rect] = imcrop(rgb);
figure;
imshow(J);
%RGB to grayscale
i = rgb2gray(J);
figure;
imshow(i);
%gaussian filter
iblur = imgaussfilt(i,0.2);
figure;
subplot(1,2,1); imshow(iblur);
subplot(1,2,2); imhist(iblur);
%histogram equalization
imh = histeq(iblur);
figure;
subplot(1,2,1); imshow(imh);
subplot(1,2,2); imhist(imh);
%edge detection
ED = edge(imh,'Canny');
figure;
imshow(ED);
%closing
se = strel('disk',4);
closeBW = imclose(ED,se);
figure;
imshow(closeBW);
%opening
se = strel('disk',3);
afterOpening = imopen(closeBW,se);
figure;
imshow(afterOpening);
%binary image
BW = imbinarize(afterOpening);
BW0=BW;
figure;
imshow(BW0);
Matt J
Matt J 2021 年 2 月 19 日 5:00
Because afterOpening is already a logical image, there is no reason here to be using imbinarize. The code you need is simply,
BW = bwpropfilt(~afterOpening,'EquivDiameter',1);
[y,x]=find(bwskel(BW));
obj=circularFit([x,y].'); %https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics
imshow(BW);
hold on
viscircles(obj.center,obj.radius);
hold off
Ahmad Najmuddin
Ahmad Najmuddin 約7時間 前
ok. thanks Matt

サインインしてコメントする。

その他の回答 (3 件)

Image Analyst
Image Analyst 2020 年 12 月 22 日
Try this:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures if you have the Image Processing Toolbox.
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 gray scale image.
folder = pwd;
baseFileName = 'crescent moon.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
% Didn't find it there. Check the search path for it.
fullFileName = baseFileName; % No path this time.
if ~exist(fullFileName, 'file')
% Still didn't find it. Alert user.
errorMessage = sprintf('Error: %s does not exist.', fullFileName);
uiwait(warndlg(errorMessage));
return;
end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
grayImage = rgb2gray(grayImage);
end
% Display the original image.
subplot(2, 3, 1);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;
%===============================================================================
% Crop the image because he supplied a screenshot with axes tick marks, etc. instead of the actual image.
grayImage = grayImage(16:649, 58:1288);
subplot(2, 3, 2);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Cropped Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
%===============================================================================
% Segment the image.
% Display the histogram
subplot(2, 3, 3);
imhist(grayImage);
grid on;
title('Cropped Image Histogram', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
% Specify a threshold. The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 60;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);
% Create a binary image
mask = grayImage > threshold;
subplot(2, 3, 4);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;
% Get rid of noise less than 4 pixels in area.
mask = bwareaopen(mask, 4);
% Now call imclose to get them to merge together.
se = true(1, 7); % Structuring element is a horizontal bar.
mask = imclose(mask, se);
% Extract the largest blob only.
mask = bwareafilt(mask, 1);
subplot(2, 3, 5);
imshow(mask);
title('Final Mask', 'FontSize', fontSize);
drawnow;
%===============================================================================
% Fit the data to a circle.
[y, x] = find(mask); % find() returns (row, column) which is (y, x)
[xCenter, yCenter, radius, a] = circfit(x, y);
%===============================================================================
% Display the original image with the circle in red.
subplot(2, 3, 6);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('(x, y) = (%.1f, %.1f). Radius = %.1f', xCenter, yCenter, radius);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Display circle in the overlay.
hold on;
viscircles([xCenter, yCenter], radius);
if xCenter >= 1 && yCenter >= 1
plot(xCenter, yCenter, 'r+', 'LineWidth', 2, 'MarkerSize', 50);
end
msgbox('Done! Thank you Image Analyst!');
%==========================================================================================
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R. A is an optional
% output describing the circle's equation:
%
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0
% by Bucher izhak 25/oct/1991
n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));
end

  5 件のコメント

表示 2 件の古いコメント
Image Analyst
Image Analyst 2020 年 12 月 28 日
I'm not Matt - I'm Image Analyst. If you have a question for Matt, reply to his message. If you meant to reply to me, say so. By the way, I don't like your code. Edge detection is almost never necessary, especially in a situation like this where you can use regular thresholding which is simpler.
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 11 日
Oh, I'm so sorry Image Analyst. I don't realize. Thank you so much for helping me :)
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 11 日
Can you adjust this coding using this original image as the input? I have try using your coding for this image but cannot get the same result. Sorry to trouble you Image Analyst.

サインインしてコメントする。


Image Analyst
Image Analyst 2021 年 1 月 12 日
Ahmad, here is the final code. It works beautifully with your supplied RGB image.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures if you have the Image Processing Toolbox.
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 gray scale image.
folder = pwd;
baseFileName = 'Rabiulawal 1439.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
% Didn't find it there. Check the search path for it.
fullFileName = baseFileName; % No path this time.
if ~exist(fullFileName, 'file')
% Still didn't find it. Alert user.
errorMessage = sprintf('Error: %s does not exist.', fullFileName);
uiwait(warndlg(errorMessage));
return;
end
end
grayImage = imread(fullFileName);
% Display the original image.
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;
% Crop the image
uiwait(msgbox('Click and drag out a box. Double click inside it to accept it.'));
[grayImage, rect] = imcrop();
hold on;
rectangle('Position', rect, 'EdgeColor', 'r', 'LineWidth', 2);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
grayImage = grayImage(:, :, 1); % Take the red channel - it has the most contrasty signal.
end
% Display the cropped gray scale image.
hFig = figure;
subplot(2, 3, 1);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;
% Remove useless stuff by commenting it out.
% % Gaussian filter & histogram equalization
% grayImage = imgaussfilt(grayImage,0.2);
% grayImage = histeq(grayImage);
% subplot(2, 3, 2);
% imshow(grayImage);
% axis on;
% impixelinfo; % Let user mouse around and see gray level.
% caption = sprintf('Postprocessing Image : %s', baseFileName);
% title(caption, 'FontSize', fontSize);
% impixelinfo;
%===============================================================================
% Segment the image.
% Display the histogram
subplot(2, 3, 2:3);
imhist(grayImage);
grid on;
title('Cropped Image Histogram', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
% Specify a threshold. The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 70;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);
% Create a binary image
mask = grayImage >= threshold;
subplot(2, 3, 4);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;
% Get rid of noise less than 4 pixels in area.
mask = bwareaopen(mask, 4);
% Now call imclose to get them to merge together.
se = true(1, 7); % Structuring element is a horizontal bar.
mask = imclose(mask, se);
% Extract the largest blob only.
mask = bwareafilt(mask, 1);
subplot(2, 3, 5);
imshow(mask);
title('Final Mask', 'FontSize', fontSize);
drawnow;
%===============================================================================
% Fit the data to a circle.
[y, x] = find(mask); % find() returns (row, column) which is (y, x)
[xCenter, yCenter, radius, a] = circfit(x, y);
%===============================================================================
% Display the original image with the circle in red.
subplot(2, 3, 6);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('(x, y) = (%.1f, %.1f). Radius = %.1f', xCenter, yCenter, radius);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Display circle in the overlay.
hold on;
viscircles([xCenter, yCenter], radius);
if xCenter >= 1 && yCenter >= 1
plot(xCenter, yCenter, 'r+', 'LineWidth', 2, 'MarkerSize', 50);
end
msgbox('Done! Thank you Image Analyst!');
%==========================================================================================
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R. A is an optional
% output describing the circle's equation:
%
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0
% by Bucher izhak 25/oct/1991
n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));
end

  3 件のコメント

Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 13 日
Thanks Image Analyst. But its only working for this image. For other crescent image it does not work. I attach another image that have thin crescent.
Image Analyst
Image Analyst 2021 年 1 月 13 日
So adapt it so it works with your other image - that's how you learn. You only gave me one image and never said there were more. If I fix this one, then you'll have another one, and another one. Anyway you said Matt's code worked because you accepted his answer so just use his code.
Since it seems very much like a one-off sort of application (i.e. you're not needing to analyze hundreds or thousands of images), why don't you just hand trace the points on these few images and then fit to a circle? See attached demo. It will be far less time consuming than trying to make a robust app that handles all kinds of situations with varying sky colors and varying cloud covers.
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 13 日
Thank you so much Image Analyst. Sorry for troubling you.

サインインしてコメントする。


Matt J
Matt J 2020 年 12 月 21 日
編集済み: Matt J 2020 年 12 月 22 日
It looks like you could just use bwareafilt to find the largest black object.

  11 件のコメント

表示 8 件の古いコメント
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 12 日
This is the coding that I add some function:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures if you have the Image Processing Toolbox.
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 gray scale image.
folder = pwd;
baseFileName = 'Rabiulawal 1439.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
% Didn't find it there. Check the search path for it.
fullFileName = baseFileName; % No path this time.
if ~exist(fullFileName, 'file')
% Still didn't find it. Alert user.
errorMessage = sprintf('Error: %s does not exist.', fullFileName);
uiwait(warndlg(errorMessage));
return;
end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
grayImage = rgb2gray(grayImage);
end
% Crop the image
[grayImage,rect] = imcrop(grayImage);
figure
imshow(grayImage);
axis on;
% Display the original image.
subplot(2, 3, 1);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;
% Gaussian filter & histogram equalization
grayImage = imgaussfilt(grayImage,0.2);
grayImage = histeq(grayImage);
subplot(2, 3, 2);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Postprocessing Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
%===============================================================================
% Segment the image.
% Display the histogram
subplot(2, 3, 3);
imhist(grayImage);
grid on;
title('Cropped Image Histogram', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
% Specify a threshold. The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 60;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);
% Create a binary image
mask = grayImage > threshold;
subplot(2, 3, 4);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;
% Get rid of noise less than 4 pixels in area.
mask = bwareaopen(mask, 4);
% Now call imclose to get them to merge together.
se = true(1, 7); % Structuring element is a horizontal bar.
mask = imclose(mask, se);
% Extract the largest blob only.
mask = bwareafilt(mask, 1);
subplot(2, 3, 5);
imshow(mask);
title('Final Mask', 'FontSize', fontSize);
drawnow;
%===============================================================================
% Fit the data to a circle.
[y, x] = find(mask); % find() returns (row, column) which is (y, x)
[xCenter, yCenter, radius, a] = circfit(x, y);
%===============================================================================
% Display the original image with the circle in red.
subplot(2, 3, 6);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('(x, y) = (%.1f, %.1f). Radius = %.1f', xCenter, yCenter, radius);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Display circle in the overlay.
hold on;
viscircles([xCenter, yCenter], radius);
if xCenter >= 1 && yCenter >= 1
plot(xCenter, yCenter, 'r+', 'LineWidth', 2, 'MarkerSize', 50);
end
msgbox('Done! Thank you Image Analyst!');
%==========================================================================================
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R. A is an optional
% output describing the circle's equation:
%
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0
% by Bucher izhak 25/oct/1991
n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));
end
Image Analyst
Image Analyst 2021 年 1 月 12 日
Yeah, you added some functions that broke it. I never told you to do a guassian blur followed by a histogram equalization. I virtually never do histogram equalization - it's almost never needed. It's like edge detection. It's something people learn about early in their image processing training and so they think it's something that should be used every time when in fact it rarely is. I'll fix the code and post it as another answer, though I thought you got it all working with the code in the Answer you accepted.
Ahmad Najmuddin
Ahmad Najmuddin 2021 年 1 月 13 日
So, the image does not need pre-processing to enhance it?

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by