histgramTest(action,varargin)

% in this test program, we calculate  the cumulative histogram in a local
%  window centered at each pixel,this local cumulative histogram can be
%  used to drive the level set for image and texture segmentation.
% Author:     Associate Prof. Yuanquan Wang,
% Affiliation: Tianjin Key Lab of Intelligent Computing and Novel Software Technology, 
%              School of Computer Science, Tianjin University of Technology, Tianjin 300191, China
% 01/20/2008
% Reference:  
%   1. Tony Chan, Selim Esedoglu, and Kangyu Ni, Histogram Based Segmentation Using Wasserstein Distances, SSVM 2007, LNCS 4485, pp. 697C708, 2007.
%   2. Kangyu Ni, Xavier Bresson, Tony Chan, Selim Esedog, Local Histogram based Segmentation using the Wasserstein Distance,
%         at: www.math.lsa.umich.edu/~esedoglu/Papers_Preprints/chan_esedoglu_ni.pdf 
%         or at :ftp://ftp.math.ucla.edu/pub/camreport/cam08-47.pdf 



function histgramTest(action,varargin)
if nargin<1,
   action='Initialize';
end;
    feval(action,varargin{:});
return;


function Initialize()
global HDmainf;

scrsz = get(0,'ScreenSize');
ww = scrsz(3);
hh = scrsz(4);
startp = ww*0.25;
endp = hh*0.25;
fw = scrsz(3)*1/2;
fh = scrsz(4)*0.5; 

HDmainf = figure('Position',[startp endp fw fh],...
    'Color', [0.9 0.9 0.9], ... 
    'NumberTitle', 'off', ...              
    'Name', 'Piecewise Constant Chan-Vese model without ReInit & with LBF- A demo', ... 
    'Units', 'pixels');

% fpath = '../images/bklimages/';   
fpath = '../images/testimages/';
orgname = 'noisyNonUniform';
fname1 = strcat(fpath,orgname);  fname = strcat(fname1,'.bmp');

dot = max(find(fname == '.'));
suffix = fname(dot+1:dot+3);
if strcmp(suffix,'pgm') | strcmp(suffix,'raw')
    e = rawread(fname);
else   e = imread(fname); end
if isrgb(e),  e = rgb2gray(e);  end
if isa(e,'double')~= 1,
    e = double(e);
end
maxe = max(max(e)');
mine = min(min(e)');
img = (e - mine)/(maxe - mine); 
Img = e;
[ysize xsize] = size(Img);

fpos = get(HDmainf,'Position');
fw = fpos(3); fh = fpos(4);
k = 1;
xpos = 0.5*(fw +50- xsize*k);
ypos = 0.5*(fh -50- ysize*k);
HDorigPic1=subplot(1,1,1);  imshow(img);
set(HDorigPic1,'Units', 'pixels','Position',[xpos ypos ysize*k xsize*k],'Units','normal');  
title('Original image');

% after read in the image, we will calculate the histogram using the
% following function , the histogram is stored in temp as a 3D array.
% one can display the local histogram at a pixel,say,(30,30) as follows:
% figure;plot(squeeze(temp(30,30,:)),'r');
winSize = [20 20];
[temp,intLength] = LocalSpectHist(Img,winSize);
disp(' for debug only ...')


function  [LSH,vectorLen] = LocalSpectHist(Img,winSize)
% first, we should know the range of the value for the whole image
minI = min(Img(:));
Img = Img - minI; % normalize the min of the image to zero
maxI = max(Img(:)); % hence, the range of the histogram will be [0, maxI]
vectorLen = round(maxI) + 1;% nonlinear transform to transform the image value to range [0 255]
if vectorLen> 256,
    Img = nonTran(Img,256);
    vectorLen = 256;
end

Img = round(Img);
pixHist(1:vectorLen) = 0;
pHist(1:vectorLen) = (1:vectorLen)-1;

%----------------------------------------------------------
[row,col] = size(Img);
LSH = zeros(row,col,vectorLen);

%----------------------------------------------------------
nr = floor(0.5*winSize(1));
nc = floor(0.5*winSize(2));
for k = 1:row,
    for l = 1:col,
        % we should extract the image window of size sinSize, and count the
        % pixel values,store in pixHist.
       
        winImg = Img(max(1,k - nr):min(row,k +nr),max(1,l - nc):min(col,l +nc));
        pixHist = histc(winImg(:),pHist);% much faster than by using loop?
        
        % four-fold loop, very slow if image size is big!!!
%         [ni,nj] = size(winImg);        
%         for ii = 1:ni,
%             for jj = 1:nj,
%                 n = winImg(ii,jj);
%                 pixHist(n+1) = pixHist(n+1) +1;
%             end
%         end 
%                
        LSH(k,l,:) = pixHist'/(ni*nj);
        pixHist(1:vectorLen) = 0;
    end
    pause(0.1);
end
% fprintf(1,'\n');

function tranImg = nonTran(Img,maxl)
tranImg = [];
theta = max(Img(:)) - min(Img(:));
temp = 1.0./( 1 + exp( -Img/theta  ) );
tranImg = (maxl - 1)*( temp - min(temp(:)))/(max(temp(:)) - min(temp(:)));