function w=OrientationBandpass(siz,oLow,oHigh) % w=orientationBandpass(size,oLow,oHigh) returns a "window", i.e. a matrix % meant to be used as a band-pass filter. The matrix size is mxn if "size" % is [m,n], and nxn if "size" is n. The matrix elements represent gain at % each freq, uniformly spaced from about -1 to 1 of Nyquist frequency (see % FREQSPACE). oLow and oHigh are the cut-off orientations in degrees (e.g. % 0 to 360). The interval will include all angles from oLow up to but not % including orientation oHigh. oHigh must be in the range [oLow,oLow+180]. % Setting oHigh=oLow+180 will produce an all-pass filter. The filter has % gain 1 in the orientation intervals [oLow,oHigh) and % [oLow+180,oHigh+180), and gain 0 outside those intervals. Note the % asymmetry: oLow is included and oHigh is excluded. Frequency zero, which % has no orientation, is always passed with gain 1. The logical OR % of filters for several contiguous bands, e.g. 0 to 20 and 20 to 40 will % equal the filter for the composite band, 0 to 40. oLow can be any finite % value, but only the range 0 to 180 is unique. o+180 is equivalent to o. % Here's a typical use, to produce orientation-bandpass noise: % noise=randn(n,n); % filter=OrientationBandpass(n,oLow,oHigh); % if any(any(filter~=1)) % skip all-pass filter % ft=filter.*fftshift(fft2(noise)); % noise=real(ifft2(ifftshift(ft))); % end % Also see Bandpass2, Bandpass, FREQSPACE. % 4/30/02 dgp wrote it, based on Bandpass2. % 5/2/02 dgp Made interval bounds asymmetric, omitting oHigh. if nargin~=3 error('Usage: w=OrientationBandpass(size,oLow,oHigh)') end if length(siz)<1 | length(siz)>2 error('"size" must be an integer or a pair of integers. Usage: w=OrientationBandpass(size,oLow,oHigh)') end if length(siz)==1 siz=[siz siz]; end if ~(isfinite(oLow) & isfinite(oHigh)) error('oLow and oHigh must be finite') end if oHighoLow+180 error(sprintf('oHigh %g must be in the range [oLow,oLow+180]. oLow is %g.',oHigh,oLow)); end % handle common special case quickly if oHigh-oLow==180 w=ones(siz); return end % call to meshgrid based on FREQSPACE.m t1 = ((0:siz(2)-1)-floor(siz(2)/2))*(2/(siz(2))); t2 = ((0:siz(1)-1)-floor(siz(1)/2))*(2/(siz(1))); [t1,t2] = meshgrid(t1,t2); a=atan2(t2,t1)*180/pi; % direction, -180 to 180 deg. t1=t1.*t1+t2.*t2; % radial frequency squared clear t2 a=mod(a,180); % orientation, 0 to 180 deg. oLow=mod(oLow,180); oHigh=mod(oHigh,180); w=ones(size(t1)); % start with all-pass filter if oHigh>=oLow d=find((a=oHigh) & t1>0); % find out-of-band frequencies else d=find((a=oHigh) & t1>0); % find out-of-band frequencies end if ~isempty(d) w(d)=zeros(size(d)); % zero the gain at those frequencies end