function [c] = savgol(z,d,k,ra,rb,theta) % function [c] = savgol(z,d,k,ra,rb,theta) % % Directional 2D Savitsky-Golay filtering with elliptical support. % The computation is done with a convolution, so the boundary of the % output will be biased. The boundary is of size floor(max(ra,rb)). % % INPUT % z Values to fit. % d Degree of fit, usually 2 or 4. % k Coefficient to return in [1,d+1], 1 for smoothing. % ra,rb Radius of elliptical neighborhood, ra=major axis. % theta Orientation of fit (i.e. of minor axis). % % OUTPUT % c Coefficient of fit. % % David R. Martin % March 2003 if d<0, error('d is invalid'); end if k<1 | k>d+1, error('k is invalid'); end ra = max(1.5,ra); rb = max(1.5,rb); ira2 = 1 / ra^2; irb2 = 1 / rb^2; wr = floor(max(ra,rb)); wd = 2*wr+1; sint = sin(theta); cost = cos(theta); % 1. compute linear filters for coefficients % (a) compute inverse of least-squares problem matrix filt = zeros(wd,wd,d+1); xx = zeros(2*d+1,1); for u = -wr:wr, for v = -wr:wr, ai = -u*sint + v*cost; % distance along major axis bi = u*cost + v*sint; % distance along minor axis if ai*ai*ira2 + bi*bi*irb2 > 1, continue; end % outside support xx = xx + cumprod([1;ai+zeros(2*d,1)]); end end A = zeros(d+1,d+1); for i = 1:d+1, A(:,i) = xx(i:i+d); end A = inv(A); % (b) solve least-squares problem for delta function at each pixel for u = -wr:wr, for v = -wr:wr, zz = zeros(wd); zz(v+wr+1,u+wr+1) = 1; yy = zeros(d+1,1); ai = -u*sint + v*cost; % distance along major axis bi = u*cost + v*sint; % distance along minor axis if ai*ai*ira2 + bi*bi*irb2 > 1, continue; end % outside support yy = cumprod([1;ai+zeros(d,1)]); filt(v+wr+1,u+wr+1,:) = A*yy; end end % 2. apply the filter to get the fit coefficient at each pixel c = conv2(z,filt(:,:,k),'same');