This code is just slightly modified version of the 2D low pass filter code of Peter Kovesi and Rob Gaddi.
The original 2D lowpass filter code and other filters are on:
http://www.owlnet.rice.edu/~elec301/Projects01/image_filt/matlab.html
Peter Kovesi’s MATLAB and Octave Functions for Computer Vision and Image Processing:
http://www.csse.uwa.edu.au/~pk/Research/MatlabFns/
Code:
% LOWPASSFILTER – Constructs a low-pass butterworth filter.
%
% usage: f = lowpassfilter(sze, cutoff, n)
%
% where: sze is a three element vector specifying the size of filter
% to construct.
% cutoff is the cutoff frequency of the filter 0 – 0.5
% n is the order of the filter, the higher n is the sharper
% the transition is. (n must be an integer >= 1).
%
% The frequency origin of the returned filter is at the corners.
%
% See also: HIGHPASSFILTER, HIGHBOOSTFILTER, BANDPASSFILTER
%
%%%%% Code History %%%%%
%
% October 1999
% Peter Kovesi pk@cs.uwa.edu.au
% Department of Computer Science & Software Engineering
% The University of Western Australia
%
% December 2001
% Modified Rob Gaddi gaddi@rice.edu
% ELEC 301
% Rice University
%
function f = lpf3(sze, cutoff, n)
if cutoff < 0 || cutoff > 0.5
error(‘cutoff frequency must be between 0 and 0.5′);
end
if rem(n,1) ~= 0 || n < 1
error(‘n must be an integer >= 1′);
end
% Modification ELEC 301 Project Group, Dec 2001
% Original code [rows, cols] = sze was not accepted by Matlab
rows = sze(1);
cols = sze(2);
hght = sze(3);
% End Alteration
% X and Y matrices with ranges normalised to +/- 0.5
% and extend them into three dimensional matrix
x = repmat( (ones(rows,1) * [1:cols] – (fix(cols/2)+1))/cols, [1 1 hght]);
y = repmat( ([1:rows]‘ * ones(1,cols) – (fix(rows/2)+1))/rows, [1 1 hght]);
% Z matrix with ranges normalized to +/0 0.5 with extention into three
% dimension.
z_pre = ([1:hght]‘ – (fix(hght/2)+1))/hght;
for j=1:hght
z(:,:,j)= z_pre(j)*ones(rows,cols);
end
radius = sqrt(x.^2 + y.^2 + z.^2); % A matrix with every pixel = radius relative to centre.
% Alteration, ELEC 301 Project Group, Dec 2001
% Original code fftshifted the filter before output. Since
% imFFT and imIFFT already shift, the output should remain low-centered.
f = 1 ./ (1.0 + (radius ./ cutoff).^(2*n)); % The filter
% End Alteration
