clear; close all; 
mex ClocalNormalizeDouble.c; % local normalization of type double
mex ClocalNormalize.c; % local normalization
%% parameters 
scaleFilter = .7;  % scale of Gabors, length = 17 pixels
numOrient = 15;  % number of orientations
saturation = 6.; % saturation level for sigmoid transformation
binSize = .2;  % binsize for computing histogram of q()
%% parameters for normalization
localOrNot = 1; % if we use local normalization or not. If not, set it to -1 
localHalfx = 20; localHalfy = 20; % the half range for local normalization, has to be quite large
thresholdFactor = .01;  % divide the response by max(average, maxAverage*thresholdFactor)
%% read in negative examples 
backgroundImageName = dir('backgroundImage/*.jpg');  % folder that contains two examples as background q(I)
backgroundNumImage = size(backgroundImageName, 1); % number of background examples
backgroundImageSize = zeros(backgroundNumImage, 2); % size of background examples
for img = 1 : backgroundNumImage
    tmpIm = imread(['backgroundImage' '/' backgroundImageName(img).name]);
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end
    backgroundI{img} = imresize(double(tmpIm), 1, 'nearest'); 
    backgroundImageSize(img, :) = size(backgroundI{img});
end
backgroundsizex = min(backgroundImageSize(:, 1)); backgroundsizey = min(backgroundImageSize(:, 2)); 
for img = 1 : backgroundNumImage
    backgroundI{img} = backgroundI{img}(1:backgroundsizex, 1:backgroundsizey); 
end   % make the sizes of the training images to be the same
%% filtering negative examples
disp(['start filtering background images']);
tic
[allFilter, allSymbol] = MakeFilter(scaleFilter, numOrient);  % generate Gabor filters 
halfFilterSize = (size(allFilter{1}, 1)-1)/2;  % half size of Gabor
allFilteredbackgroundImage = ApplyFilterfftSame(backgroundI, allFilter, localOrNot, localHalfx, localHalfy, 1, thresholdFactor);   % filter training images
disp(['filtering time: ' num2str(toc) ' seconds']);
%% compute histogram of q()
numBin = floor(saturation/binSize)+1;  % binnumbers
histog = zeros(numBin, 1);  % store F  
mex Chistogram.c;   % compile C code
disp(['start histogramming background images']);
tic
Chistogram(backgroundNumImage, numOrient, allFilteredbackgroundImage, halfFilterSize, backgroundsizex, backgroundsizey, binSize, numBin, histog, saturation);
disp(['histogramming time: ' num2str(toc) ' seconds']);
%% compute stored lambda, expectation, logZ
r = (0:(numBin-1))*binSize;
figure; plot(r, histog); 
title('background density'); xlabel('sigmoid(response)'); ylabel('density');
numStoredPoint = 50; 
storedExpectation = zeros(numStoredPoint, 1); 
Z = zeros(numStoredPoint, 1); 
for (k=1:numStoredPoint)
   lambda = (k-1.)/10.; 
   p = exp(lambda*r).*(histog'); 
   Z(k) = sum(p*binSize); p = p/Z(k); 
   storedExpectation(k) = sum(r.*p*binSize);
end
storedlambda = (0:(numStoredPoint-1))/10.; 
storedLogZ = log(Z); 
figure; plot(storedlambda, storedExpectation);
title('mean versus lambda'); xlabel('lambda'); ylabel('mean');
figure; plot(storedlambda, storedLogZ);
title('logZ versus lambda'); xlabel('lambda'); ylabel('logZ');
save 'storedExponentialModel' numStoredPoint storedlambda storedExpectation storedLogZ ...
      saturation scaleFilter numOrient halfFilterSize allFilter allSymbol ...
      localOrNot localHalfx localHalfy thresholdFactor; 
  




