Sinc filtering of an image

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',[])

[Matt J]

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]

Is convolution what he was talking about? I had no idea since I've never heard the term sinc interpolation that I can recall.

[Simon]

Hi, thanks for you answers! Why shouldn't I expect a sinc interpolation? 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), but there' no visible effect. Is it more precise now? Is there anybody, who could give me a code for filtering an image with sinc?

[Image Analyst]

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.

[Matt J]

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.

Answer 1

[Image Analyst]

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.

[Matt J]

Since convolution is happening in the frequency domain, the effect in the space domain will be a truncation/windowing of the image.

[Simon]

Sry, then it has to be a sinc filter, not an interpolation.

function SincFilter
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)+ ((2*sin(pi.*x))/(pi.*x)).*(x~=0);
B = y.*C;
figure(2); clf; hold on; 
imshow(B,'DisplayRange',[])
new_image(:,:) = ifft(B);
figure(3); clf; hold on; 
imshow(new_image,'DisplayRange',[])

Now it displays: "Error using .* Matrix dimensions must agree.

Error in SincFilter (line 34)" y is a 256x256 complex double and C a 1x41 double.... but how can I make them to have the same dimensions.

[Image Analyst]

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)

[Simon]

Hey :-)

Thanks a lot! That was really helpful for learning! My image is already gray, that's why I deleted this:

      if numberOfColorChannels > 1
        myImage = rgb2gray(myImage);

Now I tried to do this with a dicom image, but it shows me an black windows only. In another code I can display the same image without problems...

A = dicomread('example.IMA');
myImage = A(:,:);

subplot(2, 2, 1);
imshow(myImage)
title('Original Gray Scale Image', 'FontSize', fontSize)

Any idea why this doesn't work? All of the other images in this code are working.

[Image Analyst]

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]

After standing up I saw a slight shape of my Image in the black window... But why is it so dark? The results of your test:

whos A
whos myImage
min(A(:))
max(A(:))
  Name        Size              Bytes  Class     Attributes

    A         256x256            131072  uint16   

    Name           Size              Bytes  Class     Attributes

    myImage      256x256            131072  uint16              

ans =
      0
ans =
   2540

I have to convert it into a double, right?

A = dicomread('example.IMA');
myImage = A(:,:);
myImage = double(myImage) +1;

Now it is completely white...

[Image Analyst]

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]

Ah, that makes sense. Thanks for your explanation!

Answer 2

[Matt J]

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.