Showing posts with label Image compression. Show all posts
Showing posts with label Image compression. Show all posts

Friday, March 8, 2013

Wiener Deconvolution: Deriving the final formula


Recently i watched a video which talks about noise and filtering.That why i am posting this.After searching the same in wiki, i didn't see a detailed derivation of the formula. so I thought to write this and post.

Wiener deconvolution is a frequency domain technique which helps to remove the noise if we know the degradation function(H(t)) in advance. The degradation function degrades the image quality. It can be blurring like Gaussian blur or it can be motion blur.

X(x,y) -> Passes through degradation function ->  Some noise gets added here(N) -> our final signal Y(x,y)

If we multiply the Y with Wiener filter then it will provide an approximation of X. Lets call this approximation as X$ .Now am going to write down the full derivation. a part of this is available in wiki. But those who are not remembering complex number maths may find it difficult.

Mean square error between approximation and original data is = E | X(x,y) - X$(x,y) | ^2

but we know X$ is  can be obtained by G(x,y) * Y(x,y)  [ * is convolution operator ]
 = E( | X(x,y) -  G(x,y) * Y(x,y)  |^2)
But if convert all to frequency domain convolution becomes multiplication, ie a(t)*b(t) = a(f)b(f)
Also we know Y is X * H + N
= E(| X -  G ( XH + N)) |^2)
= E(| X -  GXH - GN) |^2)
= E(| (1 -  G H)X - GN) |^2)

if A and B are two complex numbers | A B |^2 = (A B)(AB)* , here * is complex conjugate.

= E(  ((1-GH)X - GN)((1-GH)X - GN)* )
= E(  ((1-GH)X - GN)((1-G*H*)X* - G*N*)) , here all * means complex conjugate

Expanding this is a painful thing. but i am not giving up. So are you ?

= E( ((1-GH)X(1-GH)* XX* - (1-GH)XG*N* - GN(1-GH)*X* + GNG*N* )
= ( (1-GH)(1-GH)* E{|X|^2} - (1-GH)G* E{XN*} - G(1-GH)* E{NX*} + GG* E{|N|^2}
wooha,
Lets assume noise and input signal has no relation,That why it is noise. So we can just ignore the terms E{XN*} ,E{NX*}.I think this is because at a particular frequency we can't find any relation between input data and noise. So the probability is null here. That why taking E( expectation) which sums the probability of values. here since we can't find any probability we can assume it 0.

We can denote S(f) = E{|X|^2} , it is the power spectrum of the complex values of the input data
We can denote N(f) = E{|N|^2} , it is the power spectrum of the complex values of the noise  data

so final equation becomes

minError(f) =  (1-GH)(1-GH)* S + GG*N
Now we need to find the minimum error which we can obtain for G. So following the basic calculus (like how we will find the min value of a function ) . careful about G* ,it is not G. So we can treat that as constant.

d(minError(f)/dG = -H(1-GH)*S + G*N = 0 ,
Now we need to find G from this. we will find it in a minute.


That is  G*N - H[1-GH]*S = 0

G*N = HS[1-GH]*
G*N = HS[1-G*H*]
G*N = S[H-G*HH*]
G*N = S[H-G*|H|^2 ]
G*N = SH-SG*|H|^2 ]
G*N = (G*) [ (1/G*)SH-S|H|^2 ]
N = [ (1/G*)SH - S|H|^2 ]
N + S|H|^2 = SH/G*
G*/SH  = 1/ (N+ S|H|^2 )
G* = SH / ( N+  S HH* ).         Take out s
G* =  SH / S( N/S+ HH* )
N/S is noise to signal ratio, We can replace it with a constant K


G* =  H / (K+ HH* )
G   = H*/ (K+H*H), H*H is |H|^2|
G   = H*/ (K+|H|^2)

Yes we reached the final answer , G = H*/(K+|H|^2).
if you know H, then you can find G, Then just multiply this with Fourier transformed input data. you will get a better picture. Thus you can remove some noise from the input data. In next post i will try to post some samples after applying G or wiener filter. Be ready for it.

Monday, December 17, 2012

Image Compression using DCT (Discrete Cosine Transform)

It has been a while , so I decided to write something from memory to keep this blog alive. Currently I am stuck with some differential equations,i really want to solve that by hand. The more i dig, it becomes more complicated.

DCT aka Discrete Cosine Transform is similar to Fourier transform but uses only Cosine functions.
Its output is just scalar values , not complex numbers as in Discrete Fourier Transform. In image processing point of view it can be used to compress data. Jpeg format uses this technique to compress data. Besides this you can use this property to any field , where you want to get only some influencing data out of your sample.

There exists different kind of kernel functions for DCT. Out of which common form is shown below.
X_k =
 \sum_{n=0}^{N-1} x_n \cos \left[\frac{\pi}{N} \left(n+\frac{1}{2}\right) k \right] \quad \quad k = 0, \dots, N-1.

I copied this from wiki. If you carefully examine you can visualize this as projecting your data on the cosine.

Let us now apply DCT over an array of numbers
result = FourierDCT[{1, 10, 2, 3, 4, 2, 5, 7, 2, 10, 10, 11, 2}]
{19.1372, -3.7222, 1.37076, 2.10949, -2.33856, 2.01418, -4.42612, 0.769566, -2.33272, -3.44945, -3.3978, 1.07119, -2.33659}

This is the result , now lets take only first 10 from result array. That is .
len = Length[result];  // assign length (13) to variable len.
compressed  = Take[result, 10]

{19.1372, -3.7222, 1.37076, 2.10949, -2.33856, 2.01418, -4.42612, 0.769566, -2.33272, -3.44945}

Now lets take Inverse DCT of this compressed data which is expected to give our valuable original data back. Note original array 'result' contained  13 numbers  So we need to fill the remaining items with 0's.

FourierDCT[PadRight[compressed, len], 3];
This will output as below.
{1.682, 8.26, 4.01, 1.549, 4.43, 2.42, 4.41, 6.876, 3.40, 7.369, 13.12, 8.47, 2.96}
Note the similarity with original data , We just did a simple compression on input data and recovered successfully. This is a lossy compression technique. Applying this on image will not reduce the overall quality of image. I mean it is still possible to understand the original image.

Let us now apply this example on an image. This is our original image, Did you able to recognize him? if you are a German you must. He is the prince of mathematicians.

This image's dimensions is {220,282} and GrayScale.

Lets now take only first 100 rows from this image and compress it


dctResult = FourierDCT[ImageData[ image, "Byte"]];
Length[dctResult];
// Here we are only taking first 100 rows from image , that is we are compressing it
compressed = Take[dctResult, {1, 100}];

Now let us recreate this compressed data again , and see how it looks like
padded = PadRight[compressed, {282, 220}];
resImage = FourierDCT[padded, 3];
Image[resImage, "byte"]



It is not that bad, still you can recognize him.
How about taking 200 rows from the original image. See the output below . It is actually good. We just saved (282-200) * 220 bytes by this simple compression. (282 is original image height, 220 is width)




I used mathematica program here. Hopes still you could get the idea.
Now again , I am asking , Did you recognize this genius ? if not he is Carl Friedrich Gauss, the Legend.