How to measure the moon crescent

Hi, anyone can help me how to measure the moon crescent using regionprops and how to remove the unwanted object in the image?

Answer 1

[Matt J]

Right now i want to measure the area of the crescent and then viscircles on it.

Using this FEX package,

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

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

[Ahmad Najmuddin]

Thanks Matt. Really helpful. Can I know how do you convert the black area of the crescent into white colour?

[Matt J]

It's just

BW=~imbinarize(Image);

[Ahmad Najmuddin]

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]

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]

ok. thanks Matt

Answer 2

[Image Analyst]

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

[Ahmad Najmuddin]

Wow! Amazing Matt. Thanks. That's what I want. Is it get the same result if I insert the original photo?

[Ahmad Najmuddin]

%upload rgb image

rgb = imread('Rabiulawal 1439.png');

[J,rect] = imcrop(rgb);

figure

imshow(J);

axis on;

%RGB to grayscale

i = rgb2gray(J);

figure;

imshow(i);

axis on;

%gaussian filter

iblur = imgaussfilt(i,0.2);

figure;

subplot(1,2,1); imshow(iblur);

subplot(1,2,2); imhist(iblur);

axis on;

%histogram equalization

imh = histeq(iblur);

figure;

subplot(1,2,1); imshow(imh);

subplot(1,2,2); imhist(imh);

axis on;

%edge detection

BW = edge(imh,'Canny');

figure;

imshow(BW);

axis on

%closing

se = strel('disk',4);

closeBW = imclose(BW,se);

figure;

imshow(closeBW);

axis on;

%opening

se = strel('disk',3);

afterOpening = imopen(closeBW,se);

figure;

imshow(afterOpening);

axis on

Please comment on my coding also Matt. I'm still new with the MATLAB so the coding so simple. I stuck at the segmentation step. The steps that I want to use consist of image acquisition, pre-processing, segmentation, feature extraction and object detection.

[Image Analyst]

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]

Oh, I'm so sorry Image Analyst. I don't realize. Thank you so much for helping me :)

[Ahmad Najmuddin]

https://www.mathworks.com/matlabcentral/answers/699060-how-to-measure-the-moon-crescent#answer_582080

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.

Answer 3

[Image Analyst]

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

[Ahmad Najmuddin]

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]

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]

Thank you so much Image Analyst. Sorry for troubling you.

Answer 4

[Matt J]

It looks like you could just use bwareafilt to find the largest black object.

https://www.mathworks.com/help/images/ref/bwareafilt.html

[Ahmad Najmuddin]

I already made the segmentation. Right now i want to measure the area of the crescent and then viscircles on it.

[Image Analyst]

No, don't use an edge detection like all those other novices. There is no justification for that. Just take your original binary image like you initially posted and use bwareafilt() like Matt suggested.

Better yet, give us your original gray scale image. Chances are you can simply just threshold it, get the points using find(), and apply the circle fitting code from the FAQ:

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F

[Matt J]

I already made the segmentation

No, it does not appear that way. A successful segmentation would mean that only the crescent is visible in your segmentation map. I see a swarm of many small black blobs that need to be removed.

[Ahmad Najmuddin]

This is my original grayscale image

[Image Analyst]

Not really. The original image would not have a white frame around it and tick marks and tick labels. Please attach the original image, not a screenshot of it after you've used imshow() or something to display it in an axes control on a figure window.

[Ahmad Najmuddin]

This is my original photo Image Analyst.

[Image Analyst]

OK, once you crop the image, it looks like my original answer should work fine.

https://www.mathworks.com/matlabcentral/answers/699060-how-to-measure-the-moon-crescent#answer_582080

Does it not?

[Ahmad Najmuddin]

I had try to add some function such as crop the image, apply gaussian filter and histogram equalization (preprocessing) to your coding but cannot get the good result. The circle not fit well with the crescent and the final image not shown the original image. Please help me Image Analyst.

[Ahmad Najmuddin]

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]

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]

So, the image does not need pre-processing to enhance it?