Salt and Pepper Image Denoising Using MATLAB

Look at the image above. If you squint, you might be able to make out what it is.

I made this image by importing a public domain image into MATLAB by typing this into the Image Processing Toolbench enabled console:

A = imread("noise.png");
A=imnoise(uint8(A),'salt & pepper',0.6);

The result is a scrambled image that doesn't convey meaningful information. Areas with random minimum and maximum value pixels characterize salt and pepper style distortion. This will either be random white and black pixels in grayscale pictures or random red, green, or blue pixels in color images. This type of problem usually happens in digital systems where a digital "1" turns into a digital "0" (or the reverse!) by mistake.

Suppose we want to write a MATLAB script to denoise our image. This post isn't about solving differential equations, so we'll leave that for another day.
Below is our Total Variation equation created by the mathematicians Rudin, Osher, and Fatemi.

This series of first and second derivatives will calculate the Minimum of the Total Variation Energy. The 2 lambda D(u-f) portion of the function acts as a fidelity factor to attempt to minimize "false positive" denoising, which would erode the original image. D is a variable for a denoising mask that we'll create by hand (well, kind of) before running our denoising function.

Let's get started.

This article assumes that you have MATLAB with the Image Processing Toolbox installed. If you're a student, you can obtain an educational license with many of the most common toolboxes preloaded for $99 at the time of writing. If not, you could use this as an informational guide for writing your algorithms in another language or try to use the Linux MATLAB clone called "Octave." I'm also assuming that you know how to operate MATLAB. If not, consider getting this book. It helped me figure out what I was doing for basic functions.

For this program, you'll manually type in the D = lines in italics into the command line to set up your mask. You'll then create a function named TV_inpaint with all the code up to the "end function" comment. You'll then run the function by typing the TV_inpaint command with the three required inputs into the command line. Good luck!
First, we will look at the code to create our denoising mask:

D = ones(size(A));
D = D - D.(A255 | A0);

This creates a mask that detects only the Salt and Pepper noises in the image. It looks for maximum or minimum values in the image and then creates a binary mask.

If you look at the mask, it'll just look like noise. That's fine. Doing this will make our image look much better after we run our script. Speaking of which...

function [u] = TV_inpaint (f, lambda, D)
% Calculate the TV energy value of an image.
% Input: Grayscale image f
% Output: TV Energy of image f
set(gcf,'color','white'); % white background
f = double(f);

This converts our input into double notation. We'll convert it back to 8 bits later.

dt = 0.5;
T = 300;
R = 255rand(size(f));
f = D.*f + (1-D).*R;
u = f;

This is our constant declaration. dt acts as our time step. Setting this higher will make the function run faster but can lead to instability. T is our maximum time to run. You can set it as high as you want and exit out of the MATLAB script with CTRL+C. R is a randomization factor that will make our mask more effective. We'll apply it to our input image f, then set our image that'll be evolved, u, equal to f.

for t = 0:dt:T
[m,n] = size(f);
% First derivatives
f_x = ( u(:,[2:n,n]) - u(:,[1,1:n-1]) ) / 2;
f_y = ( u([2:m,m],:) - u([1:m-1,1],:) ) / 2;
% Second derivative in x: f_xx
f_xx = u(:,[2:n,n]) - 2u + u(:,[1,1:n-1]);
% Second derivative in y: f_yy
f_yy = u([2:m,m],:) - 2u + u([1,1:m-1],:);
% Diagonal derivative f_xy
f_xy = (u([2:m,m],[2:n,n]) + u([1,1:m-1],[1,1:n-1])- u([1,1:m-1],[2:n,n])- u([2:m,m],[1,1:n-1])) / 4;
% Gives TV
u = u + dt.*((f_xx.f_y.^2 - 2f_x.*f_y.*f_xy + f_yy.f_x.^2)./((f_x.2 +f_y.^2).(3/2)+0.1) - 2lambda.D.(u-f));
imshow(uint8(u)); title(t);
drawnow;
end % end for
end % end function

This huge for-loop calculates the first and second derivative of our modified input image, f. It then applies the Total Variation denoising equation to our output image, u. We loop back up and do it all over again... hundreds and hundreds of times. I've written this script to do 3000 iterations. This will take a while to run but will produce a good output.
Alright, now that we've written everything, let's run our script:

TV_inpaint(A,0.2,D);

This will input A as our image under test, D as our denoising mask, and applies a Lambda factor of 0.2.

For this post, I ran this script on the image at the top and received this as my output. I think it looks pretty good, considering the source image.