function real_filt = hurlbert( gamma, beta, sigma, lambda )
% Anya Hurlbert's Regularization Approach to Color Constancy
% This function generates a filter with the indicated parameters:
%
% hurlbert( gamma, beta, sigma, lambda )
%
% gamma  (typ. 1)    - the noise reduction term
% beta   (typ. 0.01) - the weight on surface albedo changes
% sigma  (typ. 4)    - the std. deviation of gaussian constant color patch
% lambda (typ. 10)   - the weight on smooth lighting changes

% J. Watlington, 11/29/96

var = sigma^2;

M = 128;	   % Length of the filter to be generated (odd)
MM = M/2;   % Number of samples in the frequency domain

%  These are the frequencies we are sampling at
omega = ( 2 * pi * [0:M] )/M;

%  Calculate the frequency response at the sample points
p = zeros( 1, M );    % for speed
for k = 1:MM+1;
	freq = omega(k);
	p(k) = (lambda*freq^2)/(lambda*freq^2 + (1+lambda*freq^2)*(beta*exp(-(freq^2)*var) + gamma*freq^4));
end
for k = MM+2:M;
    p(k) = p(M-k+2);
end

figure;
subplot( 2, 1, 1 );
	plot( omega(1:M), p, '-' );
	axis( [ 0, pi, 0, 1.1 ] );
	xlabel( 'omega' ), ylabel( 'magnitude' );
	title( 'Frequency Response' );

text1 = [ 'gamma: ',  num2str( gamma ),  '  beta:  ', num2str( beta ) ];
text2 = [ 'lambda: ', num2str( lambda ), '  sigma: ', num2str( sigma ) ];

text( 2, 0.8, text1 );
text( 2, 0.6, text2 );

rfilt = ifft( p, M );
filt( 1:MM ) = rfilt( MM+1:M );
filt( MM+1:M ) = rfilt( 1:MM );

%  This code generates a trapezoidal window designed to leave the
%  desired filter with a good DC response.
W = 4;
window_size = M/W;
window_inc = W/M;
window = ones( 1, M );
for k = 1:window_size;
	window( k ) = k * window_inc;
	window( M - k + 1 ) = k * window_inc;
end

corr_factor = 0.23;
window( M/2 + 1 ) = 1 - corr_factor;
window( M/2 ) = 1 - 2 * corr_factor;
window( M/2 - 1 ) = 1 - corr_factor;
end

for k = 1:M;
%  In conjunction with the above code, this provides a trapezoidal window
	real_filt(k) = 1.3 * real( filt(k) ) * window(k);
%  The following line adds a Hamming window, and gooses the gain a tad
% 	real_filt(k) = 1.1 .* real(filt(k)) .* ( 0.54 - 0.46 * cos( 2*pi*(k-1)/M ));
%  The following line adds a Blackman window
%	real_filt(k) = real(filt(k)) .* ( 0.42 - 0.5 * cos( 2*pi*(k-1)/(M-1)) + 0.08 * cos( 4*pi*(k-1)/(M-1)) );
end

subplot( 2, 1, 2 );
	plot( [-MM:(MM-1)], real_filt( 1:M ), '-' );
	axis( [ -MM, MM, -0.1, 0.6 ] );
	title( 'Spatial Response' );

ffilt = fft( real_filt, 2048 );
figure;
	pomega = 0:pi/1024:pi;
	plot( pomega, abs( ffilt( 1:1025 )), '-' );
	axis( [ 0, pi/4, 0, 1.3 ] );
	xlabel( 'omega' ), ylabel( 'magnitude' );
	title( 'Freq. Response of actual filter' );