Sinc filtering of an image
Show older comments
Hi,
I want to use a sinc filter to interpolate an image. Unfortunately theres no difference between the orginal image and the new image. Could somebody help me? I don't find the mistake :-( Thanks!
function SincFilter
%read image
A = imread('example.bmp');
x1 = -20:1:20;
x = abs(x1);
org_image = A(:,:);
figure(1); clf; hold on;
imshow(org_image,'DisplayRange',[])
M = size(org_image, 1);
N = size(org_image, 2);
y(:,:) = fft(A);
% % sinc(x) = (sin(pi*x))/(pi*x);
C = (1).*(x==0)+ ((sin(pi.*x))/(pi.*x)).*(x~=0);
B = conv2(y(:,:),C, 'same');
figure(2); clf; hold on;
imshow(B,'DisplayRange',[])
new_image(:,:) = ifft(B);
figure(3); clf; hold on;
imshow(new_image,'DisplayRange',[])
5 Comments
Matt J
on 10 Dec 2012
Convolving the spectrum y with a sinc is the Fourier dual of multiplying A by a rect window. There's no reason to expect it to implement sinc interpolation.
BTW, using the expression y(:,:) instead of just y is inefficient. It forces MATLAB to create a copy of y unnecessarily.
Image Analyst
on 10 Dec 2012
Is convolution what he was talking about? I had no idea since I've never heard the term sinc interpolation that I can recall.
Image Analyst
on 10 Dec 2012
Because that's not interpolation. That's filtering. Interpolation is figuring out what values lie in between other values. You're not doing that.
I haven't run your case but it's possible that your sinc is so narrow that it's basically a delta function and so your output image would look virtually identical to your input image.
The image must be smoothed by convolving the image with the sinc function(kernel) or muliplying the fft of the image with the fft of the sinc function(kernel)
You're doing neither. You are convolving the fft of the image with a sinc, not multiplying it with one.
Also, as ImageAnalyst points out, smoothing and interpolation are different things, so your goals are not clear.
Accepted Answer
Image Analyst
on 10 Dec 2012
0 votes
It's possible that your sinc is so narrow that it's basically a delta function and so your output image would look virtually identical to your input image. Try widening your sinc.
8 Comments
Matt J
on 10 Dec 2012
Since convolution is happening in the frequency domain, the effect in the space domain will be a truncation/windowing of the image.
Image Analyst
on 10 Dec 2012
Edited: Image Analyst
on 10 Dec 2012
The sinc is the FT of the rect function. I'm still not sure what you want to do. A typical thing to do is to convolve the image with a 2D rect function in the spatial domain. Recall that convolving in the spatial domain is multiplying in the spectral domain. So convolving your image with a 2D rect function is the same as multiplying your spectrum by a sinc in the spectral domain. So if you look at the spectrum for the original image, and then look at the spectrum of the filtered image, you should see that it looks multiplied by a sinc function. HOWEVER, often the spectrum of real images has a huge DC and low frequency component so that the tails of the sinc will be affecting very small values so you might not notice. Two things to do in that case: 1) take the log of the spectral image before you display it to compress the lower frequencies and see the effect of the multiplication at the high frequencies, and 2) adjust the width of the rect function. Recall that widths operate inversely in the two domains, so that the narrower the rect function in the spatial domain, the wider the sinc in the spectral domain. And the wider the rect function in the spatial domain, the narrow the sinc in the spectral domain. You might have to play around with the width of the rect and observe what happens to the log of the image in the spectral domain. It will be a great learning experience for you. In the meantime, run this demo to see how you can use log to see the spectrum better, and to see the effect of a sinc multiplication in the Fourier domain caused by filtering with a box filter in the spatial domain:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures.
clear; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;
%2D FFT Demo
%Import images
myImage = imread('cameraman.tif');
[rows columns numberOfColorChannels] = size(myImage)
if numberOfColorChannels > 1
myImage = rgb2gray(myImage);
end
%Display images
subplot(2, 2, 1);
imshow(myImage)
title('Original Gray Scale Image', 'FontSize', fontSize)
%Perform 2D FFTs
fftA = fft2(double(myImage));
shiftedFFT = fftshift(fftA);
subplot(2, 2, 2);
imshow(real(shiftedFFT));
title('Real Part of Spectrum', 'FontSize', fontSize)
subplot(2, 2, 3);
imshow(imag(shiftedFFT));
title('Imaginary Part of Spectrum', 'FontSize', fontSize)
%Display magnitude and phase of 2D FFTs
subplot(2, 2, 4);
imshow(log(abs(shiftedFFT)),[]);
colormap gray
title('Log Magnitude of Spectrum', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Now convolve with a 2D rect function.
figure;
rectWidth = 10;
rectHeight = 5;
kernel = ones(rectHeight, rectWidth) / (rectHeight * rectWidth);
% Display it
subplot(2, 2, 1);
k = padarray(kernel, [3, 3]); % Just for display.
imshow(k, []);
axis on;
title('Kernel', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Convolve kernel (box filter) with the image
filteredImage = conv2(double(grayImage), kernel, 'same');
% Display filtered image.
subplot(2, 2, 2);
imshow(filteredImage,[]);
title('Filtered Image', 'FontSize', fontSize)
% Perform 2D FFT on the filtered image to see its spectrum.
% We expect to see a sinc multiplication effect.
% It should look like the original but with a sinc pattern overlaid on it.
fftFiltered = fft2(double(filteredImage));
shiftedFFT = fftshift(fftFiltered);
% Display magnitude of the 2D FFT of the filtered image.
subplot(2, 2, 3);
imshow(log(abs(shiftedFFT)),[]);
colormap gray
title('Log Magnitude of Spectrum - Note sinc multiplication', 'FontSize', fontSize)
Image Analyst
on 11 Dec 2012
Is it double? What do these say
whos A
whos myImage
min(A(:))
max(A(:))
Don't use semicolons so you can see what it reports to the command window.
Simon
on 11 Dec 2012
Image Analyst
on 11 Dec 2012
No, don't convert to double. The reason is that your image is 16 bit so the image display range is 0 to 65535 by default. But your brightest pixel is only 2540 - way less than 65535 - so it looks very dark. You can do this to scale it automatically so that 2540 becomes 255 instead of 65535 becoming 255:
imshow(myImage, []);
The double became all white because if the class is double, it expects it to be in the range zero to one. Zero goes to 0 and 1 goes to 255. Anything greater than 1 is treated like it's 1 - it gets clipped. Since all your pixels were greater than 1, they got clipped to the max display intensity of 255.
Simon
on 11 Dec 2012
More Answers (1)
The Fourier dual of sinc filtering by multiplication in the frequency domain is rect filtering by convolution in the space domain. If that is your goal, then it would be much simpler (2 lines of code) to work in the space domain instead,
rect1D=ones(1,n);
new_image = conv2(rect1D,rect1D,org_image, 'same');
It is also more efficient computationally to do it this way, assuming the rect window width n is not too large, which is usually the case in practice.. If you intend n to be large, you're going to destroy the image no matter which domain you work in.
Categories
Find more on Matched Filter and Ambiguity Function in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
