Introduction:
For our project we decided to tackle the task of implementing designing software to perform various image degradation and restoration. We decided that Matlab would be the best program to use, as our previous homework gave us some experience working with this software. Of course because our previous homework had covered a similar topic there had to be something that differentiated our project from the homework. We accomplish this by not using any of the built in Matlab functions for our program, but instead wrote our own functions to emulate the ones found already within the program. Starting off we first decided to implement certain noise models, which will be list later on the in the report, so we could test our filters later on in the program. After designing these types of noise we needed build some filters to clean the images up, so we went with many types of filters that our listed later on in this report also. This project while not necessarily breaking any new ground in the image processing world will help us to explore the nature what is actually happening when someone invokes a Matlab command. The hope is that someone reading this would gain a newfound respect for the work that Matlab performs in the background. Now I will list the individual modules of our program and discuss the functionality and any problems that we had while implementing them.
FREQUENCY DOMAIN FILTERING FOR IMAGE RESTORATION
BANDPASS & BANDREJECT
The process of eliminating sinusoidal noise using the frequency domain includes preprocessing, Fourier transform of the preprocessed image, creating a filter to block out unwanted frequencies, multiplying the filter and Fourier transformed image in the frequency domain, taking the inverse Fourier transform of that image, and finally to post process the image. The initial image which has some noise is defined as “f”, and the final image after post processing is complete is defined as “g”. Pre and post-processing include multiplying the each pixel/element of the images “f” or “g1” by (-1)^(x+y). After preprocessing is complete, the Fourier transform is taken and the resulting image labeled F. Now the mask can be generated by a predefined equation, or a mask that the user wishes to use. After the mask is generated, it is multiplied by F in the frequency domain. The Inverse Fourier transform is then taken of the resulting image to get ‘g1’, and finally after post-processing, the image restoration is complete. In the case of the Band-pass and Band-reject filters, the masks consist of a circular region defined by an inner and outer radius of width W. Inside the walls of the circle the signal is 1, signaling a Band-pass filter or a 0, which signals a Band-reject filter. The three types of masks that we created for both Band-pass and Band-reject were Ideal, Butterworth, and Gaussian filters. The benefits of Gaussian and Butterworth filters are that in the frequency domain, they have very nice smooth edges with little or no ripple. The Fourier transform of the Ideal Band-pass filter is a sinc function, which, has very noticeable ripples in the final image, “g”. The equations used in the program are listed below:
Ideal Band-pass/reject Filter:
If-
the distance from the origin of the frequency domain is greater than the selected frequency, Do, minus half of the width of the desired band-pass/reject filter, W and less than the selected frequency plus the width of the desired band-pass/reject filter;
Then-
the circular region defined is 1 inside the circle and zero outside for a band-pass. For a band-reject filter, inside the circular region is 0, and outside is a 1.
Gaussian and Butterworth Filters(element-wise operations)
Butterworth-
Band-reject:
Hbbr(u,v) = 1 / ( 1 + [( D(u,v)*W)/(D2(u,v) - Do2)]2N
Band-pass:
Hbbp(u,v) = 1-Hbbr(u,v)
Gaussian-
Band-reject:
Hgbr(u,v) = 1- exp(.5*(D^2(u,v)-Do^2)/D(u,v)*W))^2)
Band-pass:
Hgbp(u,v) = 1-Hgbr(u,v)
Notch Filters
Notch filters are done in the same way as the band-pass/reject filters except are filled in circles and not necessarily based at the origin. They “notch” into the surrounding mask, by becoming the opposite of the surrounding. So a notch reject, would be a small reject region in a larger pass region. This is done if there are any frequencies that are not based at the origin, but still affect the image and must be removed. One quality of the notches that is from the Fourier transform is the dual notches which are the same distance and direction from the origin, but 180 degrees rotated around the origin. We made three types of notch filters Ideal, Butterworth, and Gaussian. The equations of these functions are listed in the Matlab code section of the report.
Eliminating Atmospheric and Motion blur using Inverse Filtering
If an image is blurred there are a few ways that we took to try to fix the problem. One was to use the Gaussian Low-pass Filter, and the other was to use Inverse Filtering. Inverse filtering is the method of dividing the masked image in the frequency domain divided by the degradation function that caused the blur. This is probably not known and must be estimated. The problem with the inverse filter is that when the degradation function is really close to zero, the accuracy decreases and the resulting image is useless. A solution that we implemented is to use the Gaussian Low-pass Filter with a high radius, which seemed to give better, but not great results. The equations for the low-pass filter and the inverse filter are listed in the Matlab code.
Spatial Filtering – Eliminating Salt and Pepper Noise
Salt and pepper noise is white and black pixels which occur in random positions on a picture. To eliminate this pesky noise, there are many options. One can average each pixel with the surrounding n pixels and save that value at the original pixel. This does a good job with very little noise, but doesn’t help very much with a large amount of noise. A better solution is to take the same amount of pixels that surround each pixel, and instead of averaging them, to sort them and take the middle value. This is called Median filtering and is much better at eliminating the noise and doesn’t have a blurring effect that the averaging method does. This works very well, but there is an improvement in this method called the Adaptive Median Filter. This method checks each pixel/element, and determines if the surrounding pixels and selects the median pixel to use. Before it saves the pixel, the algorithm checks if the pixel is an impulse, and if it is increases the size of the range of pixels selected. Again the median pixel is chosen, and then checked if it is an impulse. If the median pixel is not an impulse, then the pixel is saved and the algorithm moves on to the next pixel. This algorithm gave the best results and could even make pictures that were not distinguishable from the human eye easy to see. All of these algorithms are listed in the Matlab code.
Fourier Transform:
It became evident at the beginning of this project that for some of the noise and filter functions we would need to perform a 2D fourier transform. This seemed like an easy task at first, but ended up being a somewhat difficult problem. After my first assessment of the equation, I decided to write the code using for loops and just let Matlab loose on all of the computations. Unfortunatly after letting Matlab work on a function for 15 minutes with still no result I decided to compute the number of loop iterations that Matlab was performing. This number ended up being 256^4(over 4 billion) because of the 4 for loops that were required for u, v, x, and y. Obviously this was not the way to go, so I tried to tackle the equation itself, trying to simplify it down using basic mathematical properties. This also failed, and after a certain amount of debugging I was forced to abandon this approach and go in a new direction. It was that of the Fast Fourier Transform that led to the eventual answer. Using my knowledge gained in another class(ECE 431) I was able to write up a 1D transform that would take a significantly less computational time. The only problem would be porting it over to the second dimension. I accomplished by looking at the code for my another project I have worked on, which was a verilog version of a simplified JPEG encoder/decoder. Looking at the discrete cosine transform I decided that I could perform the FFT on the horizontal part the image and then on the vertical part, the end result was a 2D fourier transform that while a little slow behave just as the Matlab command fft2.
Atmospheric Turbulence:
This noise model is based on the equation H(u,v) = e^(-k((u^2) + (v^2))(5/6), which is in the frequency domain. So using the 2D fourier transform from above I computed the transform of our input image. Then I used a for loop to compute H(u,v) and then performed a element-wise matrix multiplication on these two matrices. Of course the result was still in the frequency domain, but after taking the inverse fourier transform of the result, an atmospheric degraded image was the result.
Conclusion:
The point of this project was to program up our own versions of Matlab functions so that we could understand the commands better, and gain more understanding of the noise and filters that these commands implemented. I would say that we accomplished these goals, as our code and descriptions show that we have a very good handle on all of these concepts. Leaving this project, we definitely have a much deeper understanding of the filters and noise models, plus much more respect for Matlab and how much easier it makes image processing.