MATLAB Answers

Sajid Khan

calculating kurtosis of impulse noise

Asked by Sajid Khan
on 29 May 2013

Hi Everyone,

I want to know the kurtosis of salt and pepper noise. It's a noise which are added on the extreme values of an image due to saturation.

A paper says that it's always greater than 20, but according to my observations, it's greater then 20 when percentage of noise is less than 10. For higher percentages, it's smaller.

see the link for a histogram of impulse noise.


1 Answer

Answer by Image Analyst
on 30 May 2013

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()
		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.
		% 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.
		folder = fullfile(matlabroot, '\toolbox\images\imdemos');
		if ~isdir(folder)
		  errorMessage = sprintf('Error: The following folder does not exist:\n%s', folder);
		% 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);
		% 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);
			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);
			% 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') 
			% Let's compute and display the histogram.
			[pixelCounts grayLevels] = imhist(grayImage);
			subplot(2, 2, 4); 
			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;
			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')
				promptMessage = sprintf('%s\n\nDone with demo!', message);
		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);
	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
% "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)
		% 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);
		set(handles.txtInfo, 'String', errorMessage);
	return; % from ComputeImageMoments


Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

MATLAB Academy

New to MATLAB?

Learn MATLAB today!