A new fast denoising algorithm for auction

Here is a novel denoising algorithm. It denoises a 256x256 color image in 2.17 seconds, grayscale in 0.58 seconds (should this be plural or singular?) in Matlab.

I will give the code with explanation to the highest bidder. Minimum bid $1000 USD.

Noisy image

Denoised image

Details of the noisy image

Details of denoised image:
Notice the texture of the wig is almost perfectly maintained.

Nostaligia

It is snowing in Toronto and it is depressing and gloomy to look outside my window. This is a perfect time to think about the past. When I was a child we used to have a film camera. Tricky contraption, if you ask me. The light would be imparted on a film, we would call it a 'negative'. This film is then developed into a photograph. For some reason, I was more fascinated by the negative. Somehow I would think that all my negativity is captured in that. So, on this gloomy day I found an old picture of mine and made its negative. Can you figure out a one line code to do this?

Me and my alter ego

A quick divergence to the cutest Laplacian in the world.

Let $\mathbf{u}=(u_1, u_2)$ be a vector valued function. Then the divergence operator: $\mbox{div }\mathbf{u}:=\frac{du_1}{dx}+\frac{du_2}{dx}$ is a very useful thing to have on our side. Again, I will leave this little code upto you. It is really a one liner!

Note, that if $\nabla u=(u_x, u_y)$ then $\mbox{div}(\nabla u)$ gives the Laplacian $\Delta u$ of the image $u$. This could be useful.

Here is the Laplacian of Boo displayed between -50 to 50.
i.e. use the commands to display the Laplacian d :
>> m=-50; M=50; figure; imshow((d-m)/(M-m));

Laplacian of Boo (-50 to 50)

Here is the original.

Boo, the cutest dog in the world

The edge of darkness is upon us!

I assume that now we know how to find edges in a given image u. It's simple, just take the derivative!

Today, I want to find the smooth regions. Well, it's kind of silly, as now that we have edges, the rest of the image is kind of edgeless. Right... but I would like to assign a number from 0 to 1 to it ... where 1 indicates a flat region and zero indicates an edge, in other words I want to make the edges appear dark.

There are many ways to do it. One of the ways to do it is to look at the function $$g(x, y):=\frac{1}{\sqrt{1+|K*\nabla u(x, y)|}}.$$Where, K is your favourite smoothing kernel. (I will leave it to you to code this as an exercise. It is a four liner.)

Staring at this function is a refreshing activity, something that I love to do in my spare time.

Boo in the edge of darkness

Lenna in the edge of darkness

Barbara in the edge of darkness

Boo, the smoothest dog in the world!

Yesterday night I was too sleepy to code, so I wrote a code for filtering with a kernel that is more than 3x3. Btw, I still like my kernel to be of odd sized and square.

Let's use the image of Boo, the cutest dog in the world and make it the smoothest dog in the world. By using a 5x5 averaging kernel.


>> h=ones(5)/25;
>> u=double(imread('boo.jpg'));
>> uf=ifilter(u, h);m=min(min(min(uf))); M=max(max(max(uf))); figure; imshow((uf-m)/(M-m));

Boo, the cutest dog in the world!

Boo, after smoothing with a 5x5 kernel becomes the smoothest dog in the world!

Reflections on convolutions and animal testing

Convolution is a very tool important in image processing. There are inbuilt code to perform convolution in Matlab. But while solving many problems in image processing we assume Neumann boundary conditions. For this reason we reflected the boundary with our cute little code expand, while we performed our algorithms in the interior of the image. Finally we got rid of the extra image boundary by our shrink program. Many a times filtering is realized by performing a convolution of some kernel K with the image u.

In a discrete case, the convolution K*u is essentially the "masking process" as described in my blog about the Sobel filtering.

For example, the Laplacian operator, or the Sobel operator we have discussed before could be coded elegantly using Matlab's inbuilt filter2 command (which uses conv2).

I don't like two things about using using these two in-built commands for filtering.

1. Matlab uses "zero padding" in these in built commands, whereas we need a reflecting boundary. i.e. it adds zeros on the boundary of an image.
2. These commands can be used for only matrices.

Well, these are really trivial issues, and we can take care of them by writing our own little function say ifilter to avoid these issues.

This code performs a convolution of the image u with any 3x3 kernel h with Neumann boundary conditions. i.e. This code alone can find gradients, perform Sobel, Prewitt, Roberts filtering, or any other filtering with 3x3 window in a very nice, elegant way ! 

Of course you can write a more generalized code that is not restricted to 3x3 window. I won't do it because it is 3.00 am and I can hardly stay awake.

Let's see which filter should I chose... Let's design a filter that picks up derivatives at 45 degrees angle.

          0     1     0

h =   -1     0     1

          0    -1     0

We can normalized the above filter by dividing it by 2.

Let us do some animal testing now. I want to use an image of a Boston terrier.

Boston terrier

Here is the result of the filter, after using the following commands:

>> u=imread('bostonterrier.png');

>> h=(1/2)*[0 -1 0; -1 0 1; 0 1 0];

>> uf=ifilter(u, h);m=min(min(min(uf))); M=max(max(max(uf))); imshow((uf-m)/(M-m));

Directional derivative

               
Laplacian dog

Last time we used Laplacian, lets use ifilter to find the Laplacian of the Boson terrier.


>> h=[0 1 0; 1 -4 1; 0 1 0];
>> uf=ifilter(u, h);m=min(min(min(uf))); M=max(max(max(uf))); imshow((uf-m)/(M-m));

Laplacian dog

Prewitt operator

Here is the classical Prewitt operator implemented using ifilter. 

>> hx=(1/3)*[-1 -1 -1; 0 0 0; 1 1 1]


   -0.3333   -0.3333   -0.3333
         0         0         0
    0.3333    0.3333    0.3333


>> hy=(1/3)*[-1 0 1; -1 0 1; -1 0 1]


   -0.3333         0    0.3333
   -0.3333         0    0.3333
   -0.3333         0    0.3333


>> ux=ifilter(u, hx); uy=ifilter(u, hy); prewitt=sqrt(ux.^2+uy.^2);
>> m=min(min(min(uf))); M=max(max(max(prewitt))); imshow((prewitt-m)/(M-m));

Prewitt dog

Heat equation removes noise ... and many other things

Heat equation $\frac{\partial u}{\partial t}=\Delta u$ is extremely useful for smoothing images,
where $u(x, 0)=f(x)$ i.e. $u$ is initialized to the given original image $f$, the image to be smoothed, we also impose Nuemann boundary conditions. Smoothing an image removes noise in image $f$. The problem with heat equation is though that it can not distinguish between noise and useful features like edges. Plus you never know when to stop! It just keeps on smoothing.

The boundary condition that we could use is Neumann boundary condition. i.e. the normal derivatives at the boundary is zero. The way to maintain this is to reflect the boundary, run the heat equation in the interior of the image, then shrink it back to the original size and repeat.

The Laplacian operator $\Delta u=u_{xx}+u_{yy}$ is easy to discretize:
$\Delta u \equiv Lu=(u_{i, j-1}+u_{i-1, j}+u_{i, j+1}+u_{i+1, j}-4u_{i, j})/h^2.$

It's a five point operator,

$$
Lu=\frac{1}{h^2}\left[\begin{array}{ccc}0 & 1 & 0 \\1 & -4\,\,\,\,\,\, & 1 \\0 & 1 & 0\end{array}\right]
$$
Before we perform the heat equation on an image $f$, we should see how the image really looks like.


>> f=imread('eyes.jpg');
>> figure; imshow(f);



This is what we get...

Image f

Let us try to see view the image as a graph of a function. This image being an RGB image we have to change it to grayscale before viewing it.

>> figure; mesh(double(rgb2gray(f)));

This command gives the following output.

Image f viewed as a graph of a function

In the above image, the magnitude of the grayscale is given by with blue indicating small magnitude and red color indicating large magnitide. This is a very "spiky" image. If we run heat equation on this image we should get a smoother image.

Here is the code I wrote for the heat equation. By now I assume that you have been convinced that the for loops are bad, so one can write the code heat_2fast without for loops.

If you write the following in the command window:

>> u=heat_2fast(f, 1, 0.25, 1); figure; imshow(u/255);

Output of the heat equation at t=1.

Let us try to see view the image as a graph of a function:

>> v=rgb2gray(u/255);

>> mesh(v);

Result of heat equation, it makes the image smooth !

The reason why heat equation is called a diffusion equation

Simply put heat equation diffuses the function f. As a matter of fact it diffuses the function f  iso-tropically. i.e. equally in all directions.  Let us see what happens if we run the heat equation for t =1000. This is what I got as a result:

The result of heat equation after a long time

This really flattened things out, huh? If we see this as a graph of a function we get the following.

Heat equation FTW!

This is what happens if you run the heat equation for t=1000. If you keep it running for more time eventually things become really flat like Kansas.