Try my demo:
% Get the mean gray level, standard deviation, skew, and kurtosis from the histogram bin values.
% Processes every image in the Image Processing demos folder.
function ComputeImageMomentsDemo()
try
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures.
clear; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
fontSize = 20;
% Change the current folder to the folder of this m-file.
if(~isdeployed)
cd(fileparts(which(mfilename)));
end
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
% User does not have the toolbox installed.
message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
if strcmpi(reply, 'No')
% User said No, so exit.
return;
end
end
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
if ~isdir(folder)
errorMessage = sprintf('Error: The following folder does not exist:\n%s', folder);
uiwait(warndlg(errorMessage));
return;
end
% Check for bmp, tif, jpg, and png images.
filePattern = fullfile(folder, '*.tif');
imageFiles1 = dir(filePattern);
filePattern = fullfile(folder, '*.png');
imageFiles2 = dir(filePattern);
filePattern = fullfile(folder, '*.bmp');
imageFiles3 = dir(filePattern);
filePattern = fullfile(folder, '*.jpg');
imageFiles4 = dir(filePattern);
imageFiles = [imageFiles1; imageFiles2; imageFiles3; imageFiles4];
numberOfImageFiles = length(imageFiles);
% Bail out if there are no images in the folder.
if numberOfImageFiles == 0
errorMessage = sprintf('Error: The following folder does not contain any tif images:\n%s', folder);
uiwait(warndlg(errorMessage));
return;
end
% Preallocate space for a moments matrix.
moments = zeros(numberOfImageFiles, 4);
for k = 1 : numberOfImageFiles
% Get the baseFileName and extension.
baseFileName = imageFiles(k).name;
% Get the extension separately.
[f b extension] = fileparts(baseFileName);
% Get rid fo the dot and convert to upper case.
extension = upper(extension(2:end));
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
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);
% Let's compute moments only for gray scale images.
% If we read in a color image, convert to grayscale by taking the green channel.
if numberOfColorBands >= 2
grayImage = grayImage(:,:,2);
end
% Display the original gray scale image.
subplot(2, 2, 3);
imshow(grayImage, []);
axis on; % Show pixel coordinates along the axis (edges of the image).
caption = sprintf('Gray Scale %s-Format Image', extension);
title(caption, 'FontSize', fontSize);
if k == 1
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.
% Give a name to the title bar.
set(gcf,'name','Image Moments Demo','numbertitle','off')
end
% Let's compute and display the histogram.
[pixelCounts grayLevels] = imhist(grayImage);
subplot(2, 2, 4);
bar(pixelCounts);
grid on;
title('Histogram of Gray Scale Image', 'FontSize', fontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.
% Compute the image moments.
[meanGL stdDev skew kurtosis] = ComputeImageMoments(grayLevels, pixelCounts);
% Plot the image moments for all images up through this one.
moments(k, 1) = meanGL;
moments(k, 2) = stdDev;
moments(k, 3) = skew;
moments(k, 4) = kurtosis;
subplot(2,1,1);
plot(moments, 'LineWidth', 3);
title('Intensity Moments of Gray Scale Image', 'FontSize', fontSize);
grid on;
legend('Mean', 'Standard Deviation', 'Skewness', 'Kurtosis');
message = sprintf('There are %d pixels in %s (#%d of %d).\nThe mean gray level is %.2f.\nThe standard deviation of the gray levels is %.2f.\nThe skewness of the gray levels is %.2f.\nThe kurtosis of the gray levels is %.2f.',...
sum(pixelCounts), baseFileName, k, numberOfImageFiles, meanGL, stdDev, skew, kurtosis);
if k < numberOfImageFiles
promptMessage = sprintf('%s\n\nDo you want to Continue processing,\nor Cancel to abort processing?', message);
button = questdlg(promptMessage, 'Continue', 'Continue', 'Cancel', 'Continue');
if strcmp(button, 'Cancel')
break;
end
else
promptMessage = sprintf('%s\n\nDone with demo!', message);
msgbox(promptMessage);
end
end % Read in a standard MATLAB gray scale demo image.
catch ME
errorMessage = sprintf('Error in ComputeImageMoments().\nThe error reported by MATLAB is:\n\n%s', ME.message);
fprintf(1, '%s\n', errorMessage);
uiwait(warndlg(errorMessage));
end
return; % from ComputeImageMomentsDemo
%------------------------------------------------------------------------------------------------------
% Computes first, second, third, and fourth central moments of the gray levels.
% Get the mean gray level, standard deviation, skew, and kurtosis from the histogram bin values.
% Note: gray level moments are different than spatial moments which are more like rotational moments of intertia.
% Uses formulas from http://itl.nist.gov/div898/handbook/eda/section3/eda35b.htm
% "Skewness is a measure of symmetry, or more precisely, the lack of symmetry.
% A distribution, or data set, is symmetric if it looks the same to the left and right of the center point.
% The skewness for a normal distribution is zero, and any symmetric data should have a skewness near zero.
% Negative values for the skewness indicate data that are skewed left and positive values for the skewness
% indicate data that are skewed right. By skewed left, we mean that the left tail is long relative to the right tail.
%
% Kurtosis is a measure of whether the data are peaked or flat relative to a normal distribution.
% That is, data sets with high kurtosis tend to have a distinct peak near the mean,
% decline rather rapidly, and have heavy tails. Data sets with low kurtosis tend
% to have a flat top near the mean rather than a sharp peak. A uniform distribution would be the extreme case."
function [meanGL stdDev skew kurtosis] = ComputeImageMoments(GLs, pixelCounts)
try
% Get the number of pixels in the histogram.
numberOfPixels = sum(pixelCounts);
% Get the mean gray lavel.
meanGL = sum(GLs .* pixelCounts) / numberOfPixels;
% Get the variance, which is the second central moment.
varianceGL = sum((GLs - meanGL) .^ 2 .* pixelCounts) / (numberOfPixels-1);
% Get the standard deviation.
stdDev = sqrt(varianceGL);
% Get the skew.
skew = sum((GLs - meanGL) .^ 3 .* pixelCounts) / ((numberOfPixels - 1) * stdDev^3);
% Get the kurtosis.
kurtosis = sum((GLs - meanGL) .^ 4 .* pixelCounts) / ((numberOfPixels - 1) * stdDev^4);
catch ME
errorMessage = sprintf('Error in ComputeImageMoments().\nThe error reported by MATLAB is:\n\n%s', ME.message);
uiwait(warndlg(errorMessage));
set(handles.txtInfo, 'String', errorMessage);
end
return; % from ComputeImageMoments
