Notch filter to remove a small area around the fundamental frequency and the first few harmonics

I have an image

After taking the log of one plus the magnitude of the transform of the this image to observe the bright spots due to the stripes. I got this

TO make it easier I used the thermal version of the image and got this

What I have to do is :-

a. Design a “notch” filter to remove a small area around the fundamental frequency and the first few harmonics.(or in other words is to remove the 'dark red strips on the coke bottle'). To remove the white spots.

b. Notch filter must be symmetric about the origin of the frequency plane.

I thought of removing the values of these pixels above a particular threshold but it is not giving me exactly what I wanted.

So this is the code I have

    f=imread('Coke.tif'); % Reading the image
    [M,N]=size(f); % getting the dimensions of the image
    F=fft2(double(f)); % taking the fast fourier transform of the image
    u=0:(M-1); %Setting the coordinates X axis
    v=0:(N-1); % Setting the coordinates Y axis
    idx=find(u>M/2); % looking for X axis origin
    u(idx)=u(idx)-M; % setting the X axis origin
    idy=find(v>N/2); % Looking fo Y axis origin
    v(idy)=v(idy)-N; % Setting the X axis origin
    [V,U]=meshgrid(v,u); % shifting the FFT of the image so that the origin is at the center of the image 
    D=sqrt(U.^2+V.^2); % Taking the square root
    NF = ones(ir,ic); % setting initial coordinates to be zero
    NF(21:27,45:63) = 0; % coordinates for the bright spots
    NF(256-27:256-21,256-63:257-45) = 0;% coordinates for the bright spots
    NF(47:53,63:81) = 0; % coordinates for the bright spots
    NF(256-53:256-47,256-81:257-63) = 0;% coordinates for the bright spots
    NF(73:79,81:99) = 0; % coordinates for the bright spots
    NF(256-79:256-73,256-99:257-81) = 0;% coordinates for the bright spots
    NF(99:105,99:117) = 0; % coordinates for the bright spots
    NF(256-105:256-99,256-117:256-99) = 0;% coordinates for the bright spots
    G = F.*NF; % Apply the notch filters to the Fourier spectrum of the image
    imshow(uint8(log(1+abs(G))),[]);
    g=real(real(ifft2(double(G)))); % getting the inverse fourier transform
    imshow(f),figure,imshow(g,[ ]); % display the original image and the fourier transform of the image

but this is giving me an error

Error using .* Matrix dimensions must agree.

Error in Notchfilter (line 59) G = F.*NF; % Apply the notch filters to the Fourier spectrum of the image

Any help please?

Answer 1

[Alejandro Arrizabalaga]

Your code is not OK. ir and ic are not defined and you seem not to do iffshift, even though you say something in your comments. I think that you forgot to paste some part. By the way, I do not get the error you say. This is what I did:

   [M,N]=size(f); % getting the dimensions of the image
    F=fft2(double(f)); % taking the fast fourier transform of the image
    u=0:(M-1); %Setting the coordinates X axis
    v=0:(N-1); % Setting the coordinates Y axis
    idx=find(u>M/2); % looking for X axis origin
    u(idx)=u(idx)-M; % setting the X axis origin
    idy=find(v>N/2); % Looking fo Y axis origin
    v(idy)=v(idy)-N; % Setting the X axis origin
    [V,U]=meshgrid(v,u); 
    % shifting the FFT of the image so that the origin is at the center of the image 
    F=fftshift(F);
    D=sqrt(U.^2+V.^2); % Taking the square root
    NF = ones(size(F,1),size(F,2)); % setting initial coordinates to be zero
    NF(21:27,45:63) = 0; % coordinates for the bright spots
    NF(256-27:256-21,256-63:257-45) = 0;% coordinates for the bright spots
    NF(47:53,63:81) = 0; % coordinates for the bright spots
    NF(256-53:256-47,256-81:257-63) = 0;% coordinates for the bright spots
    NF(73:79,81:99) = 0; % coordinates for the bright spots
    NF(256-79:256-73,256-99:257-81) = 0;% coordinates for the bright spots
    NF(99:105,99:117) = 0; % coordinates for the bright spots
    NF(256-105:256-99,256-117:256-99) = 0;% coordinates for the bright spots
    G = F.*NF; % Apply the notch filters to the Fourier spectrum of the image
    G=ifftshift(G);
    figure;imshow(uint8(log(1+abs(G))),[]);
    g=real(ifft2(double(G))); % getting the inverse fourier transform
    figure;imshow(f),figure,imshow(g,[ ]); % display the original image and the fourier transform of the image