function [allfiltered] = applyfilterfftsameLocalNormPerOri(I, allfilter)
% applyfilterfftsameLocalNormPerOri Filter image by a bank of filters, with per-orientation local normalization.  
% 
% I: input images
% allfilter: filter bank
% threshold: filter response is set to zero if below threshold

nimage = size(I, 2);    % number of training images
norient = size(allfilter, 2);     % number of orientations
h = (size(allfilter{1}, 1)-1)/2;  % half size of filters
[sx, sy] = size(I{1});  % size of images
allfiltered = cell(nimage, norient);  % filtered images

fftI = cell(1,nimage);
for i = 1 : nimage
   fftI{i} = fft2(I{i}, sx+h+h, sy+h+h);
end
fftf = cell(1,norient);
for o = 1 : norient
   fftf{o} = fft2(allfilter{1, o}, sx+h+h, sy+h+h); 
end

mask = ones(sx,sy);
mask([1:h sx-h:sx],:) = 0; mask(:,[1:h sy-h:sy]) = 0;

for i = 1 : nimage
   fftIi = fftI{i}; 
   for o = 1 : norient
       fftfo = fftf{o}; 
       out = ifft2(fftIi.*fftfo); 
       filtered = out(h+1:h+sx, h+1:h+sy); 
       re = real(filtered); im = imag(filtered); 
       energy = re.*re + im.*im;  % compute the local energy
       ave = mean(mean( energy(h+1:end-h,h+1:end-h) ));
       allfiltered{i, o} = energy / ave;
   end
end


radius = 20; % 20 for good results
aveFilter = ones(2*radius+1,2*radius+1)/((2*radius+1)^2);
for i = 1:nimage
    for o = 1:norient
        energy = allfiltered{i,o};
        energy([1:2 end-1:end],:) = 0; energy(:,[1:2 end-1:end]) = 0;
        localMeanMap = imfilter(energy,aveFilter,'replicate');
        allfiltered{i,o} = sigmoidOnOriginalResponse(energy) .* energy .* mask ./ (localMeanMap + 1e-10);
    end
end

%% sigmoid transform
function v = sigmoidOnOriginalResponse(v)
v = 1 ./ ( 1 + exp( 1 * -( v - 3 ) ) );
