function [contrast_detection_threshold, adaptation_luminance] = local_adapt_detect( L, L_surround, peak_sensitivity, ppd )
%LOCAL_ADAPT_DETECT
%
% local_adapt_detect(L) computes the map of approximate contrast detection
% thresholds and the map of local adaptation luminance for the given input
% luminance image L in cd/m^2.
%
% local_adapt_detect(L, L_surround) uses L_surround as the surround
% luminance outside the image L. Default: the mean of the 4 corners of L:
% mean( [L(1,1) L(1,end) L(end,1) L(end,end)] ).
%
% local_adapt_detect(..., peak_sensitivity) sets the peak sensitivity of
% the TVI function in photopic conditions. The peak sensitivity is defined
% as the inverse of the detection threshold (L/DeltaL). For example,
% peak_sensitivity=100 means that the observer is able to detect a
% difference of 1 cd/m^2 (DeltaL) at 100 cd/m^2 (L). Default: 6, which
% corresponds to the detection threshold (DeltaL/L) of approx. 17%.
%
% local_adapt_detect(..., ppd) sets the pixels per degree of the input
% image L. Default: 45
%
% [contrast_detection_threshold] = local_adapt_detect(...) returns the
% adaptation map of contrast detection thresholds given as Weber fractions:
% DeltaL/L. For example, when the map value is 0.01 and the luminane of the
% input pixel is 100 cd/m^2, the observer can detect a difference of
% 0.01 * 100 cd/m^2 = 1 cd/m^2. 
%
% Please note that model is not meant to be a universal detection model and
% the thresholds need to be adjusted for each application with the
% "peak_sensitivity" parameter. If you need a universal detection model,
% check HDR-VDP (http://hdrvdp.sourceforge.net/).
%
% [contrast_detection_threshold, adaptation_luminance] =
% local_adapt_detect(...) additionally returns the adaptation luminance map
% in cd/m^2.
%
% Note that the model is valid only in the photopic levels of luminance, 
% above 1 cd/m^2. 
%
% The full description of the local adaptation and contrast detection model
% can be found in:
%
% Vangorp, P., Myszkowski, K., Graf, E.W., and Mantiuk, R.K.
% A Model of Local Adaptation
% ACM Transactions on Graphics 34, 6 (Proc. ACM SIGGRAPH Asia 2015)
%
% If you use this code in your research, please cite the publication above.
%
% Project web page: http://resources.mpi-inf.mpg.de/LocalAdaptation/
%
% This code us not optimized for speed. Please do NOT use it for
% benchmarking.
%
% Tested on: Matlab R2009b with the Image Processing Toolbox
%            Octave 3.8.2 with the image package
%
% Authors: Peter Vangorp <peter.vangorp@gmail.com>
%          Rafal K. Mantiuk <mantiuk@gmail.com>
%
% License: MIT

%% The MIT License (MIT)
% 
% Copyright (c) 2015 Peter Vangorp and Rafal K. Mantiuk
% 
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
% 
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

%% The code
% Default options
if ~exist('peak_sensitivity', 'var') || isempty(peak_sensitivity)
    peak_sensitivity = 6;
end
if ~exist('L_surround', 'var') || isempty(L_surround)
    L_surround = mean( [L(1,1) L(1,end) L(end,1) L(end,end)] );
end
if ~exist('ppd', 'var') || isempty(ppd)
    ppd = 45;
end

% Best model parameters
p_nl1 = [0.030594263300886 1.798038555035897 0.072142302635550 0.010000000000000]; % non-linearity 1
sigma_nr_1 = 0.428078089769146; % Gaussian 1
p_nl2 = [0.030205776178761 2.934114110750746 0.210878838078038 0.084607721580785]; % non-linearity 2
sigma_nr_2 = 0.082417606884968; % Gaussian 2
alpha = 0.653547088004559; % summation weight
nr_exp = 1; % Naka-Rushton exponent
pupil_size = 4; % in mm

% Model
L_otf = filterEyeOTF( L, ppd, L_surround, pupil_size ); % Glare / inter-ocular scatter function

adaptation_luminance = alpha*inv_nonlin(fast_gauss(fw_nonlin(L_otf,p_nl1),sigma_nr_1*ppd),p_nl1) ...
    + (1-alpha)*inv_nonlin(fast_gauss(fw_nonlin(L_otf,p_nl2),sigma_nr_2*ppd),p_nl2);

sens_fix = 0.0332; % This is the default lowest detection threshold

contrast_detection_threshold = tvi_csf_hdrvdp( min(max(L_otf, 1e-5), 1e8) )/sens_fix/peak_sensitivity ./ L;

te = naka_rushton_te( L_otf, adaptation_luminance, nr_exp ); % threshold elevation due to mal-adaptation
contrast_detection_threshold = contrast_detection_threshold .* te;

end

%% Helper functions
function L1 = filterEyeOTF( L, ppd, L_surround, pupil_size )
rho2 = create_cycdeg_image( [size(L,1) size(L,2)]*2, ppd );
mtf_filter = exp( -( rho2 / (20.9 - 2.1*pupil_size)).^(1.3-0.07*pupil_size) );
L1 = fast_conv_fft( L, mtf_filter, L_surround );
end

function D = create_cycdeg_image( im_size, pix_per_deg )
nyquist_freq = 0.5 * pix_per_deg;
half_size = floor(im_size/2);
odd = mod( im_size, 2 );
freq_step = nyquist_freq./half_size;
if( odd(2) )
    xx = [ linspace( 0, nyquist_freq, half_size(2)+1 ) linspace( -nyquist_freq, -freq_step(2), half_size(2) ) ];
else
    xx = [ linspace( 0, nyquist_freq-freq_step(2), half_size(2) ) linspace( -nyquist_freq, -freq_step(2), half_size(2) ) ];
end
if( odd(1) )
    yy = [ linspace( 0, nyquist_freq, half_size(1)+1 ) linspace( -nyquist_freq, -freq_step(1), half_size(1) ) ];
else
    yy = [ linspace( 0, nyquist_freq-freq_step(1), half_size(1) ) linspace( -nyquist_freq, -freq_step(1), half_size(1) ) ];
end
[XX YY] = meshgrid( xx, yy );
D = sqrt( XX.^2 + YY.^2 );
end

function Y = fast_conv_fft( X, fH, pad_value )
pad_size = (size(fH)-size(X));
fX = fft2( padarray( X, pad_size, pad_value, 'post' ) );
Yl = real(ifft2( fX.*fH, size(fX,1), size(fX,2) ));
Y = Yl(1:size(X,1),1:size(X,2));
end

function Y = fast_gauss( X, sigma )
ks = [size(X,1) size(X,2)]*2;
NF = pi;
xx = [linspace( 0, NF, ks(2)/2 ) linspace(-NF, -NF/(ks(2)/2), ks(2)/2) ];
yy = [linspace( 0, NF, ks(1)/2 ) linspace(-NF, -NF/(ks(1)/2), ks(1)/2) ];
[XX YY] = meshgrid( xx, yy );
K = exp( -0.5*(XX.^2 + YY.^2)*sigma^2 );
Y = zeros( size(X) );
for cc=1:size(X,3)
    Y(:,:,cc) = fast_conv_fft( X(:,:,cc), K, 'replicate' );
end        
end

function L_out = fw_nonlin( L_in, p )
nodes = 4;
l_lut = linspace( -1, log10(5000), nodes+1 );
delta = l_lut(2)-l_lut(1);
v_lut = cumsum( [0 p*delta] );
l = min(max( log10(L_in), l_lut(1) ), l_lut(end) );
L_out = interp1( l_lut, v_lut, l, 'cubic' );
end

function L_out = inv_nonlin( L_in, p )
nodes = 4;
l_lut = linspace( -1, log10(5000), nodes+1 );
delta = l_lut(2)-l_lut(1);
v_lut = cumsum( [0 p*delta] );
l_lut_dense = linspace( l_lut(1), l_lut(end) );
v_lut_dense = interp1( l_lut, v_lut, l_lut_dense, 'cubic' );
L_out = 10.^interp1( v_lut_dense, l_lut_dense, min(max( L_in, v_lut_dense(1) ), v_lut_dense(end) ), 'linear' );
end

function T = tvi_csf_hdrvdp( L_adapt )
p5 = 30.162;
p6 = 4.0627;
p7 = 1.6596;
p8 = 0.2712;
S = p5 * ( (p6./L_adapt).^p7 + 1).^(-p8);
T = L_adapt./S;
end

function Te = naka_rushton_te( L, La, n )
Te = (L.^n + La.^n).^2./(4*(L.*La).^n); 
end
