function [f] = oeFilter(sigma,support,theta,deriv,hil,vis) % function [f] = oeFilter(sigma,support,theta,deriv,hil,vis) % % Compute unit L1-norm 2D filter. % The filter is a Gaussian in the x direction. % The filter is a Gaussian derivative with optional Hilbert % transform in the y direction. % The filter is zero-meaned if deriv>0. % % INPUTS % sigma Scalar, or 2-element vector of [sigmaX sigmaY]. % [support] Make filter +/- this many sigma. % [theta] Orientation of x axis, in radians. % [deriv] Degree of y derivative, one of {0,1,2}. % [hil] Do Hilbert transform in y direction? % [vis] Visualization for debugging? % % OUTPUTS % f Square filter. % % David R. Martin % March 2003 nargchk(1,6,nargin); if nargin<2, support=3; end if nargin<3, theta=0; end if nargin<4, deriv=0; end if nargin<5, hil=0; end if nargin<6, vis=0; end if numel(sigma)==1, sigma = [sigma sigma]; end if deriv<0 | deriv>2, error('deriv must be in [0,2]'); end % Calculate filter size; make sure it's odd. hsz = max(ceil(support*sigma)); sz = 2*hsz + 1; % Sampling limits. maxsamples = 1000; % Max samples in each dimension. maxrate = 10; % Maximum sampling rate. frate = 10; % Over-sampling rate for function evaluation. % Cacluate sampling rate and number of samples. rate = min(maxrate,max(1,floor(maxsamples/sz))); samples = sz*rate; % The 2D samping grid. r = floor(sz/2) + 0.5 * (1 - 1/rate); dom = linspace(-r,r,samples); [sx,sy] = meshgrid(dom,dom); % Bin membership for 2D grid points. mx = round(sx); my = round(sy); membership = (mx+hsz+1) + (my+hsz)*sz; % Rotate the 2D sampling grid by theta. su = sx*sin(theta) + sy*cos(theta); sv = sx*cos(theta) - sy*sin(theta); if vis, figure(1); clf; hold on; plot(sx,sy,'.'); plot(mx,my,'o'); %plot([sx(:) mx(:)]',[sy(:) my(:)]','k-'); plot(su,sv,'x'); axis equal; ginput(1); end % Evaluate the function separably on a finer grid. R = r*sqrt(2)*1.01; % radius of domain, enlarged by >sqrt(2) fsamples = ceil(R*rate*frate); % number of samples fsamples = fsamples + mod(fsamples+1,2); % must be odd fdom = linspace(-R,R,fsamples); % domain for function evaluation gap = 2*R/(fsamples-1); % distance between samples % The function is a Gaussian in the x direction... fx = exp(-fdom.^2/(2*sigma(1)^2)); % .. and a Gaussian derivative in the y direction... fy = exp(-fdom.^2/(2*sigma(2)^2)); switch deriv, case 1, fy = fy .* (-fdom/(sigma(2)^2)); case 2, fy = fy .* (fdom.^2/(sigma(2)^2) - 1); end % ...with an optional Hilbert transform. if hil, fy = imag(hilbert(fy)); end % Evaluate the function with NN interpolation. xi = round(su/gap) + floor(fsamples/2) + 1; yi = round(sv/gap) + floor(fsamples/2) + 1; f = fx(xi) .* fy(yi); % Accumulate the samples into each bin. f = isum(f,membership,sz*sz); f = reshape(f,sz,sz); % zero mean if deriv>0, f = f - mean(f(:)); end % unit L1-norm sumf = sum(abs(f(:))); if sumf>0, f = f / sumf; end