How can I implement a low-pass Butterworth filter in Matlab?

Question

From this answer, I know how to create a High-pass Butterworth filter.

From this video, I know that, lowpasskernel = 1 - highpasskernel.

So, I created the following Low-pass Butterworth Filter,

function [out, kernel] = butterworth_lp(I, Dl, n)
    height = size(I, 1);
    width  = size(I, 2);

    [u, v] = meshgrid(-floor(width/2):floor(width/2)-1,-floor(height/2):floor(height/2)-1);

    % lp_kernel = 1 - hp_kernel
    kernel = 1 - butter_hp_kernel(u, v, Dl, n);

    % fft the image
    I_fft_shifted = fftshift(fft2(double(I)));

    % apply lowpass filter
    I_fft_shift_filtered = I_fft_shifted .* kernel;

    % inverse FFT, get real components
    out = real(ifft2(ifftshift(I_fft_shift_filtered)));

    % normalize and cast
    out = (out - min(out(:))) / (max(out(:)) - min(out(:)));
    out = uint8(255*out);

function k = butter_hp_kernel(u, v, Dh, n)
    uv =  u.^2+v.^2;
    D = sqrt(uv);
    frac = Dh./D;
    p = 2*n;
    denom =  frac.^p;
    k = 1./denom;

Output

This isn't a low-pass filter output.

So, what is the issue with my code?

Answer 1

Okay. I have solved the problem by following the following formula (Page #8/48),

Output

Source Code

butter_lp_kernel.m

function f = butter_lp_f(u, v, Dl, n)
    uv = u.^2+v.^2;
    Duv = sqrt(uv);
    frac = Duv./Dl;
    denom = frac.^(2*n);
    f = 1./(1.+denom);

function k = butter_lp_kernel(I, Dl, n) 
    Height = size(I,1); 
    Width = size(I,2); 

    [u, v] = meshgrid( ...
                    -floor(Width/2) :floor(Width-1)/2, ...
                    -floor(Height/2): floor(Height-1)/2 ...
                 ); 

    k = butter_lp_f(u, v, Dl, n);

ifftshow.m

function out = ifftshow(f)
    f1 = abs(f);
    fm = max(f1(:));
    out = f1/fm;
end

butterworth_lpf.m

function [out1, out2] = butterworth_lpf(I, Dl, n)

    Kernel = butter_lp_kernel(I, Dl, n);

    I_ffted_shifted = fftshift(fft2(I));

    I_ffted_shifted_filtered = I_ffted_shifted.*Kernel;

    out1 = ifftshow(ifft2(I_ffted_shifted_filtered));

    out2 = ifft2(ifftshift(I_ffted_shifted_filtered));
end

Answer 2

You didn't correctly copy the formula in the high pass filter:

 uv =  u.^2+v.^2;
 D = uv.^0.5;
 frac = Dh./D;
 p = 2*n;
 denom =  frac.^p;
 k = 1./(1+denom); --> you forgot the "1+"

New output: