% Demo to add "salt and pepper" noise to a gray scale image, then restore the image by removing this noise
% with a modified median filter that acts only on the noise pixels
% and not upon "good" non-noise pixels.
clc;	% Clear command window.
clear;	% Delete all variables.
close all;	% Close all figure windows except those created by imtool.
imtool close all;	% Close all figure windows created by imtool.
workspace;	% Make sure the workspace panel is showing.
fontSize = 15;

% Read in a standard MATLAB demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'coins.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);
% Display the original image.
subplot(2, 3, 1);
imshow(grayImage);
title('Original Gray Scale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize')); 

% Generate a noisy image with salt and pepper noise.
noisyImage = imnoise(grayImage,'salt & pepper', 0.05);
subplot(2, 3, 2);
imshow(noisyImage);
title('Image with Salt and Pepper Noise', 'FontSize', fontSize);

% Median Filter the image:
medianFilteredImage = medfilt2(noisyImage, [3 3]);
% Display the image.
subplot(2, 3, 3);
imshow(medianFilteredImage);
title('Median Filtered Image', 'FontSize', fontSize);

%----------------------------------------------------------------------------------------------------------
% Two ways to identify the noise: (1) assume it's above or below some limits - like it's just pure 0 or pure 255
% or (2) identify it as some number of gra levels away from the median value in a neighborhood.
% Method %2 is good if the noise is not pure salt and pepper noise but has values other than 0 and 255.
% Use either method #1 OR method #2, not both.  So comment one of them out!!!
%----------------------------------------------------------
% Method #1: darker or brighter than some specified levels.
darkLevel = 0;
brightLevel = 255;
% Find the noise.  It will have a gray level of either 0 or 255.
noiseMap = (noisyImage <= darkLevel | noisyImage >= brightLevel);
%----------------------------------------------------------
% Method #2: darker or brighter than some specified levels from the median filtered image.
% differenceThreshold = 64; % Whatever you want.
% Find the noise.  It will have a gray level more than differenceThreshold away from the median value in a neighborhood.
% noiseMap = imabsdiff(grayImage, noisyImage) >= differenceThreshold;
%----------------------------------------------------------

% Display the image.
subplot(2, 3, 4);
imshow(noiseMap);
title('Noise Map', 'FontSize', fontSize);

% Get rid of the noise by replacing with median.
noiseFreeImage = noisyImage; % Initialize
noiseFreeImage(noiseMap) = medianFilteredImage(noiseMap); % Replace.
% Display the image.
subplot(2, 3, 5:6);
imshow(noiseFreeImage);
title('Restored Image', 'FontSize', fontSize);
