pass - python remove noise from signal
Removing sinusoidal noise with Butterworth filter (2)
I'm trying to remove the sinusoidal noise in this image:
Here is its DFT spectrum (after applying log and arbitrary intensity scaling):
I already have a Butterworth filter to apply to this image. It will knock out the mid-frequency peaks. I'm taking care to scale it from [0..255] to [0..1.0] after loading. Here's the filter:
The results aren't great:
- Why is there still a significant amount of noise remaining in the image?
- Why is the result darker than the original image? The filter clearly does not touch the DC term, so I'd expect the average intensity to be the same.
- Why is the filter only taking out some of the peaks? It comes from a textbook so I'm inclined to believe that it is correct, but there are other peaks in the spectrum -- are they also part of the noise? I tried removing them using concentric filters, but it didn't do much good and darkened the image beyond recognition.
I've taken the image (cropped) and filter from the book Digital Image Processing by Gonzalez and Woods. In their example, the periodic noise is completely removed by the filtering, and the mean intensity of the image remains the same.
My source code for loading the image and filter, DFT, filtering, IDFT is below:
import cv def unshift_crop(comp, width, height): result = cv.CreateImage((width, height), cv.IPL_DEPTH_8U, 1) for x in range(height): for y in range(width): real, _, _, _ = cv.Get2D(comp, x, y) real = int(real) * ((-1)**(x+y)) cv.Set2D(result, x, y, cv.Scalar(real)) return result def load_filter(fname): loaded = cv.LoadImage(fname, cv.CV_LOAD_IMAGE_GRAYSCALE) flt = cv.CreateImage(cv.GetSize(loaded), cv.IPL_DEPTH_32F, 2) width, height = cv.GetSize(loaded) for i in range(width*height): px, _, _, _ = cv.Get1D(loaded, i) #cv.Set1D(flt, i, cv.Scalar(px/255.0, 0)) cv.Set1D(flt, i, cv.Scalar(px/255.0, px/255.0)) return flt if __name__ == '__main__': import sys fname, filt_name, ofname = sys.argv[1:] img = cv.LoadImage(fname, cv.CV_LOAD_IMAGE_GRAYSCALE) width, height = cv.GetSize(img) src = cv.CreateImage((width*2, height*2), cv.IPL_DEPTH_32F, 2) dst = cv.CreateImage((width*2, height*2), cv.IPL_DEPTH_32F, 2) cv.SetZero(src) for x in range(height): for y in range(width): px, _, _, _ = cv.Get2D(img, x, y) px = float(px) * ((-1) ** (x+y)) cv.Set2D(src, x, y, cv.Scalar(px, 0)) cv.DFT(src, dst, cv.CV_DXT_FORWARD) flt = load_filter(filt_name) cv.Mul(dst, flt, src) cv.DFT(src, dst, cv.CV_DXT_INV_SCALE) result = unshift_crop(dst, width, height) cv.SaveImage(ofname, result)
There was a bug in the original source where the filter imaginary components were loaded as zero. This is what caused the result image to appear darker than it really was. I've fixed this and commented the related line.
Using the fixed source and the filter that @0x69 provided (yeah, I know it's not really a Butterworth filter, but at this stage I'm happy to try anything), this is the result:
Better than what I had to start with, but still not as good as I'm hoping for. Can anyone beat this? I suspect putting more notches in to take out the remaining peaks may be of some benefit.
I've contacted the author. This is their response:
The problem is that the image used in the experiment was floating point, while the one shown in the book (and the original provided in the downloads) are 8 bits. This is required for printing, etc.
In order to duplicate the experiment, you have to start with the noise-free image and then add your own noise to it.
I remember playing around with this image during a course in image processing a couple of years ago, and I got the same results as you.
I don't know how the authors of the textbook got the image they display in the book, but they must have done something more then apply that Butterworth filter. As you mention, there are more peaks, so it is possible that they applied more Butterworth filters to remove these.
However, the mean of the image did stay the same for me. Have you tried calculating the mean of the two images and comparing? It could be that it is just the scaling when displaying that causes the darker image.
I think you make it bit Complex. Just Do it on Matlab if you like.
It gives you pretty good results.
% Question: Filtration in Frequency Domain im = imread('applo_noisy.tif'); FT = fft2(double(im)); FT1 = fftshift(FT);%finding spectrum %imtool(abs(FT1),); m = size(im,1); n = size(im,2); t = 0:pi/15:2*pi; xc=(m+150)/2; % point around which we filter image yc=(n-150)/2; r=200; r1 = 40; xcc = r*cos(t)+xc; ycc = r*sin(t)+yc; xcc1 = r1*cos(t)+xc; ycc1 = r1*sin(t)+yc; mask = poly2mask(double(xcc),double(ycc), m,n); % Convert region-of-interest polygon to mask mask1 = poly2mask(double(xcc1),double(ycc1), m,n); % generating mask for filtering % a=51; % b=5; % a(b) = 0; % 51 0 0 0 0 mask(mask1)=0; FT2=FT1; FT2(mask)=0;%cropping area or bandreject filtering output = ifft2(ifftshift(FT2)); imtool(output,);