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.

Sunday, February 24, 2013

Histogram Equalization : A simple technique for improving low contrast images

I have been busy. Not getting much time to give enough attention to this. I am really worried about this and this makes me unhappy. But life must goes on. I can do what i want but i am not able to do.. It is an interesting thought. I am sure many of us still have this.(just had a wine)

What is histogram , it is a frequency showing some information in image. Due to some factors based on the equipment or lighting or whatever it is possible sometimes taken image will contain colors that uses near by gray values. That is is the image's gray levels are not really stretched to required level.

See some images here
These are images taken on low light or artificially generated. if you examine this these uses intensities that are around near by spectrum.That is why our eye is not able to perceive the amount of detail we want to interpret.  Histogram equalization is a very simple technique to stretch out this contrast by using a cumulative distributive function.

This CDF function is very simple, it means for the pixel with highest value it assigns a probability of 1, because it is cumulative  That is for the brightest pixel value( say 255) it sums probabilities of all pixels from 0 to 255 and it must be 1 it will be (for other pixel values it sum ups probabilities up to that pixel's value). So our CDF function is a monotonically increasing function.

I am attaching the mathematica code and results. so i don't have to explain.

yimage = Import["f:\\imtests\\lowContrast1.jpg"];
grayImg = ColorConvert[myimage, "GrayScale"];
mydata = ImageData[grayImg, "byte"];
temp = ImageDimensions[grayImg];
width = temp[[1]];
height = temp[[2]];
numberOfPixels = width*height;
histo = Array[0 &, 256];

(* compute histogram *)
m = mydata[[2]][[138]];
For[ i = 1 , i <= height, i++,
  For[ j = 1 , j <= width, j++,
    m = mydata[[i]][[j]];
    m = m + 1;
    histo[[m]] ++;
    ];
  ];

(* create cummlative distributive function to map pixels *)
cdf = Array[0 &, 256];
For[ i = 1, i <= 256, i++,
  sum = 0.0;
  For[ j = 1, j <= i, j++,
 
   sum = sum + histo[[j]] / numberOfPixels;
   ];
  sum = sum * 255;
  cdf[[i]] = sum;
  ];

(* plot cdf *)
ListLinePlot[cdf]

(* equalize histogram *)

orignalData = mydata;
(* change pixels in original image based on cdf function *)
For[ i = 1 , i <= height, i++,
  For[ j = 1 , j <= width, j++,
    m = mydata[[i]][[j]] + 1;
    mydata[[i]][[j]] = cdf[[m]]
    ];
  ];
Image[orignalData, "byte"]
Image[mydata, "byte"]


Outputs

CDF  function




Original image






Output image


See how good is the equalized image. but do not expect Histogram Equalization to perform this in all conditions. In general It is good for x-ray images. One other problem is it will also amplify the noise. So it better to remove the noise before passing to this.

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.