4mar: Denoising

One of the problems of performing a low pass filter to remove noise, is that edges get blurred. There are multiple advanced techniques perform denoising of images while preserving the sharpness of edges. In this activity I try to develop simple techniques to perform edge preserving denoising, based on methods taught in class.

Activity 1a: Edge preserving denoising using variable filters

If you refer to last's week work, then you will find how we can use vertical and horizontal Sobel filters to generate edges and then combine them to get magnitude of an edge at each point. Now we want to filter our image with a low pass Gaussian filter. But the strength of the filter depends on its spread or variance. We could preserve edges, if we could make the filter weaker around edges, but stronger in flat regions. So we want the variance of the filter to be a function of the edge magnitude map that we generate. In general we want variance of impulse response to be small (ie variance of frequency response is large (all pass like character)), when the pixel has high magnitude edge and vice versa. In other words, if we are at an edge pixel, apply minimal filtering else apply strong filtering.

Below is the main loop of the code. I am just showing the code that represents the green part of the block diagram

%generate thresholds

t1 = quantile(mag(:), 0.9);

t2 = quantile(mag(:), 0.7);

t3 = quantile(mag(:), 0.3);

for i = 1:sz(1)

for j = 1:sz(2)

filter = fspecial('gaussian', r, selectvariance(mag, i, j, t1, t2, t3)); %generate gaussan filter

%do pointwise multiplication of image patch and filter window

for x = -(r-1)/2 : (r-1)/2

for y = -(r-1)/2 : (r-1)/2

if (i + y) >= 1 && (i + y) <= sz(1) && (j + x) >= 1 && (j + x) <= sz(2) %check bounds

result(i, j) = result(i, j) + filter(y + ((r-1)/2) + 1, x + ((r-1)/2) + 1) * noisyimg(i+y, j+x);

end

end

end

end

end

%a heuristic function that returns an appropriate variance based on the magnitude of edge

function var = selectvariance(mag, i, j, t1, t2, t3)

var = 1;

if (mag(i,j) > t1)

var = 0.2; return;

end

if (mag(i,j) <= t1 && mag(i,j) > t2)

var = 0.6; return;

end

if (mag(i,j) <= t2 && mag(i,j) > t3)

var = 0.8; return;

end

if (mag(i,j) <= t3)

var = 1; return;

end

end

For testing I choose a fractal image fractal formations have lots of details and edges in many different scales.

Note that in the simple filter case, I used 5x5 filter window with variance = 1. In edge preserving filter case I used variances 0.2, 0.6, 0.8 and 1 depending on edge strength.

Note that the edge preserving filter removes noise as good as the simple filter in flat regions (eg the background), but it is more sharper than the simple filtered image, as it does less filtering near the edges. This aspect is clearly seen in the detailed small triangles of the fractal.

Activity 1b: Edge preserving denoising by ignoring edges when filtering

Another simple way to preserve edges when filtering is to first do LPF an then Canny to find edges. Then we start our filtering operation. If the pixel has more than 2 pixels as edge pixels in its neighbourhood, we do not perform the filtering on that point, else we proceed with filtering. This is related to the technique described above, but is a bit faster.

for i=1:sz(1)

for j = 1:sz(2)

[n, edgepoints] = neighbourhood(im3, i, j, r); %find the number of edge points in the neighbourhood

if (n > 2)

result(i, j) = noisyimg(i, j); %many edge points around this pixel, so don't perform filtering

else

result(i, j) = multiplypatch(noisyimg, filter2, i, j, 5); %perform filtering

end

end

end

%auxiliary function to perform filtering

function val = multiplypatch(noisy, filter, i, j, r) %assume r (window length) is odd

val = 0;

sz = size(noisy);

for x = -(r-1)/2 : (r-1)/2

for y = -(r-1)/2 : (r-1)/2

if (i + y) >= 1 && (i + y) <= sz(1) && (j + x) >= 1 && (j + x) <= sz(2)

val = val + noisy(i+y, j+x) * filter(y + ((r-1)/2) + 1, x + ((r-1)/2) + 1);

end

end

end

end

%auxiliary function to find number of pixels that are edge points in the neighborhood of a pixel

function [n, edgepoint] = neighbourhood(im, i, j, r) %assume r (window length) is odd

n = 0;

if (mod(r,2) == 0)

r = r+1;

end

edgepoint = [];

sz = size(im);

for x = -(r-1)/2 : (r-1)/2

for y = -(r-1)/2 : (r-1)/2

if (i + y) >= 1 && (i + y) <= sz(1) && (j + x) >= 1 && (j + x) <= sz(2)

if (im(i+y, j+x) == 1)

n = n + 1;

edgepoint = [edgepoint x y];

end

end

end

end

end

The modified technique gives better results as the denoised image looks crisper. Also the edge detection algorithm working on this denoised image captures more edges than the simple LPF denoised image. The difference in the details lost in simple filtering is highlighted in the edge images.

Activity 2:

In this activity I describe what I did for the homework assignment. The objective is to enhance the images that suffer from specific problems by applying appropriate tools and make the result visually pleasing.

Image 1

Below is the image and its histogram

It is clear from the histogram that the image does not have enough dynamic range. The pixels are mostly concentrated towards the bright side, hence the image appears washed out. So we can rectify that with a histogram equalization. Below is the function for histogram equalizaion given an image and the original histogram.

%inputs: image and histogram of image

function eq = histeq(img, h)

sz = size(img);

eq = zeros(sz);

prob = h/(sz(1)*sz(2)); %probability table

cdftable = zeros(1,256);

cdftable(1) = prob(1);

for i = 2:256

cdftable(i) = prob(i) + cdftable(i-1);

end

%generate equalised image

for i = 1:sz(1);

for j = 1:sz(2)

eq(i,j) = uint8(floor(cdftable(img(i,j)) * 255.0));

end

end

eq = eq/max(max(eq));

h = imhist(eq); %plot new histogram

figure;plot(h);

figure; plot(cdftable);

end

We get the following result from histogram

Another very simple, but effective way to enhance this image would be to normalize it, that is map the min value to 0 and the max value to 1 and everything in between is mapped by a linear map. Below is the function and the results

function normimg = normalize(img)

normimg = img - min(min(img)); %gets rid of negative values

normimg = normimg/max(max(normimg)); %scales in the range 0 to 1

end

We could also use log transform for this image

The edge detected on the histogram equalised image is better than the edges detected on the original image. It is especially visible in the edges of the clouds, which is not detected in the original image, but is detected for the equalised image. Canny detector is used.

The flow diagram for this enhancement is as follows

Image 2

The next image looks quite dark.

Applying histogram equalization we get the following image:

The image still looks noisy and it looks like the noise is white noise. Therefore a Gaussian filter might be of help here. Below are the results from applying a Gaussian LPF of size 11x11 and variance 7

Now Let us compare Canny edge detectors results on the original image, the equalised image and then the final filtered image.

As we can see above, the original image is not easy to detect edges on because of its small dynamic range. The equalised image, shows lots of edges, but due to the texture of thegrass, lot of the edges are probably spurious. So filtering helps to reduce local perturbations, and the edge detector performs better on the filtered image and its detected edge looks less noisy.

The flow diagram for this enhancement is:

Image 3

The next image:

It is a bit tricky because at first glance the image looks corrupted with salt and pepper noise, so a simple median filtering will suffice. Below is the image processed with median filter and the edge detected through Canny.

Since Salt and pepper noise sets some pixels to black (0) and some to white (255), it might be a good idea to scan the image and apply median filtering only to those points which have 0 or 255 as their pixel values. This would save some time as we don't have to apply the median to every pixel, and this is also a good heuristic as salt and pepper noise only sets pixels to 0 and 255, so other pixels should not undergo filtering ideally. So to achieve that I used this function

function result = saltpepperfilter(img, window)

sz = size(img);

radius = (window-1)/2

for i = 1:sz(1)

for j = 1:sz(2)

%only if pixel is 0 or 255 apply median filter

if (img(i,j) == 255 || img(i,j) == 0)

%this if else structure makes sure the window is within bounds of image size

if (i - radius < 1)

ilow = 1;

else

ilow = i - radius;

end

if (j - radius < 1)

jlow = 1;

else

jlow = j - radius;

end

if (i + radius > sz(1))

ihigh = sz(1);

else

ihigh = i + radius;

end

if (j + radius > sz(2))

jhigh = sz(2);

else

jhigh = j + radius;

end

subset = img(ilow:ihigh, jlow:jhigh); %select a window around the pixel

subsetvector = subset(:);

result(i,j) = median(subsetvector);

else

result(i,j) = img(i,j);

end

end

end

end

But on applying this kind of filtering, after the salt and pepper noise was removed, it turned out that the image was also corrupted with white noise. The presence of white noise was not evident in the original picture. Also applying median filtering wholesale is not a good idea (due to the presence of white noise). Therefore after selectively applying median filtering to remove salt and pepper noise, I applied Gaussian LPF to remove the white noise.

Let us now compare the edges detected by applying Canny on simple median filtered image vs (selective median filtering + Gaussian LPF) image. As marked in yellow, some edges are lost in simple median filtering. This is because the image is corrupted with two types of noise and we should apply correct filters in correct order to get rid of the noise one by one. Note that applying selective salt and pepper filter first is a good idea as salt and pepper noise is easy to identify and remove, unlike white noise. Also if we had applied the Gaussian filter first without removing the salt and pepper noise, then the filter would have smeared the salt and pepper noise.

The summarizing flow diagram for this enhancement is:

Image 4

Here is the next image.

It is white noise corrupted so I apply Gaussian LPF. But in this image, the woman's dress and the table cloth have a lot of high frequency content (stripes). LPF blurs the stripes. Hence I filter the low pass image with Difference of Gaussian filter to get edge information and add it back to the image to enhance the edges.

Maybe the images do not show too much difference, but if you click on the image and expand it, the stripes on the the pants and scarf can be seen to have been enhanced in the second picture. Thus the second picture has a crisper look as the edge information was partially restored. The steps are summarised below.

Image 5

The next image is a colour image corrupted in all 3 channels with independent salt and pepper noise. The 3 channels are shown individually too.

After applying median filtering to each of the channels and recombining them we get the following results

We still have some noise left over, so let us apply the median filter again for a few more times

After a few rounds almost all the noise is removed. Now let us apply edge detection on individual channels as well as after converting the image to grayscale. We observe that almost all he edges are similar in the 3 channels and also in the grayscale image

The summary for this enhancement is shown below.