function DocBinarize
% user interactive: user chooses a part of the image where Otsu threshold
% provides training image
% it reads image, calls recursive subdivision function, then performs TPS
% interpolation, and finally generates binary image

close all;
clear all;

global X Y T hmin wmin vec_gt int_step rec_cnt;
X=[]; Y=[]; T=[]; hmin=16; wmin=hmin; int_step=5;

warning off;

% read image
ORG_filename = 'H02.bmp';


% gray scale image to be binarized
im = imread(ORG_filename);
if size(im,3) ==3                            % Only get a grey scale image
    im = im(:,:,1);
end
figure(1), imagesc(im),colormap(gray),axis image;
title('Create a rectangle enclosing some text: Click on the top-left and bottom-right corners');
[x,y] = ginput(2);
x=round(x);y=round(y);
figure(1),hold on, plot([x(1),x(2),x(2),x(1),x(1)],[y(1),y(1),y(2),y(2),y(1)],'r','LineWidth',2),hold off; drawnow;
% x=[74; 238]; y=[51; 118];
% x=[85;311]; y=[420; 493];
im_sml = im(y(1):y(2),x(1):x(2));
th = graythresh(double(im_sml)/255);
im_gt = im_sml>th*255;
[vec_gt, binloc_gt] = Training_GetWdthHist(im_gt); % Get the histogram

im=double(im);%/255;
% call recursive subdivision program
dist_min=1e10;
Imin = min(im(:));
Imax = max(im(:));
th_min=255;
disp('Estimating initial threshold value...');
for th=Imin+10:int_step:Imax-10,
    vec_sub = Training_GetWdthHist(im>=th);
    dist_vec = Testing_HistogramDist(vec_sub, vec_gt);
    if dist_min>dist_vec,
        dist_min=dist_vec;
        th_min = th;
    end
    disp(['Threhold value: ' num2str(th)]);
end
disp('Initial estimation done.');
figure(1),hold on;
rec_cnt=0;
RecursiveThreshold(im,1,1,th_min,dist_min);
plot(X,Y,'gs','MarkerFaceColor','g','MarkerSize',10);
hold off;
% call thin plate spline
disp('Doing spline interpolation...');
[h,w]=size(im);
TPS  = Testing_ThinPlateSpline([Y X], h, w, T);
figure(2), imagesc(im>TPS),colormap(gray),axis image;
disp('done.');


function RecursiveThreshold(I,x0,y0,th_min,dist_min)
% I: input image
global X Y T hmin wmin vec_gt int_step rec_cnt;
rec_cnt=rec_cnt+1;
disp(['Number of calls to recursive subdivision: ' num2str(rec_cnt)]);
wid=5;
[h,w] = size(I);
if h>w,
    hby2 = round(h/2);
    if hby2>=hmin && w>=wmin,
        left_dist_min=1e10;
        Imin = min(min(I(1:hby2,:)));
        Imax = max(max(I(1:hby2,:)));
        left_th_min=255;
        for th=Imin+10:int_step:Imax-10,
            vec_sub = Training_GetWdthHist(I(1:hby2,:)>=th);
            dist_vec = Testing_HistogramDist(vec_sub, vec_gt);
            if left_dist_min>dist_vec,
                left_dist_min=dist_vec;
                left_th_min = th;
            end
        end
        if left_dist_min<=dist_min,
            RecursiveThreshold(I(1:hby2,:),x0,y0,left_th_min,left_dist_min);
        else
            T=[T; th_min; th_min];
            X=[X; x0+w/4; x0+3*w/4];
            Y=[Y; y0+h/4; y0+h/4];
            rectplot([x0 x0+w-1],[y0 y0+hby2-1],wid);
        end
        right_dist_min=1e10;
        Imin = min(min(I(hby2:end,:)));
        Imax = max(max(I(hby2:end,:)));
        right_th_min=255;
        for th=Imin+10:int_step:Imax-10,
            vec_sub = Training_GetWdthHist(I(hby2:end,:)>=th);
            dist_vec = Testing_HistogramDist(vec_sub, vec_gt);
            if right_dist_min>dist_vec,
                right_dist_min=dist_vec;
                right_th_min = th;
            end
        end
        if right_dist_min<=dist_min,
            RecursiveThreshold(I(hby2:end,:),x0,y0+hby2-1,right_th_min,right_dist_min);
        else
            T=[T; th_min; th_min];
            X=[X; x0+3*w/4; x0+w/4];
            Y=[Y; y0+3*h/4; y0+3*h/4];
            rectplot([x0 x0+w-1],[y0+hby2-1 y0+h-1],wid);
        end
    else
        T=[T; th_min; th_min; th_min; th_min];
        X=[X; x0+w/4; x0+3*w/4; x0+3*w/4; x0+w/4];
        Y=[Y; y0+h/4; y0+h/4; y0+3*h/4; y0+3*h/4];
        rectplot([x0 x0+w-1],[y0 y0+h-1],wid);
    end
else
    wby2 = round(w/2);
    if wby2>=wmin && h>=hmin,
        left_dist_min=1e10;
        Imin = min(min(I(:,1:wby2)));
        Imax = max(max(I(:,1:wby2)));
        left_th_min=255;
        for th=Imin+10:int_step:Imax-10,
            vec_sub = Training_GetWdthHist(I(:,1:wby2)>=th);
            dist_vec = Testing_HistogramDist(vec_sub, vec_gt);
            if left_dist_min>dist_vec,
                left_dist_min=dist_vec;
                left_th_min = th;
            end
        end
        if left_dist_min<=dist_min,
            RecursiveThreshold(I(:,1:wby2),x0,y0,left_th_min,left_dist_min);
        else
            T=[T; th_min; th_min];
            X=[X; x0+w/4; x0+w/4];
            Y=[Y; y0+h/4; y0+3*h/4];
            rectplot([x0 x0+wby2],[y0 y0+h-1],wid);
        end
        right_dist_min=1e10;
        Imin = min(min(I(:,wby2:end)));
        Imax = max(max(I(:,wby2:end)));
        right_th_min=255;
        for th=Imin+10:int_step:Imax-10,
            vec_sub = Training_GetWdthHist(I(:,wby2:end)>=th);
            dist_vec = Testing_HistogramDist(vec_sub, vec_gt);
            if right_dist_min>dist_vec,
                right_dist_min=dist_vec;
                right_th_min = th;
            end
        end
        if right_dist_min<=dist_min,
            RecursiveThreshold(I(:,wby2:end),x0+wby2-1,y0,right_th_min,right_dist_min);
        else
            T=[T; th_min; th_min];
            X=[X; x0+3*w/4; x0+3*w/4];
            Y=[Y; y0+h/4; y0+3*h/4];
            rectplot([x0+wby2-1 x0+w-1],[y0 y0+h-1],wid);
        end
    else
        T=[T; th_min; th_min; th_min; th_min];
        X=[X; x0+w/4; x0+3*w/4; x0+3*w/4; x0+w/4];
        Y=[Y; y0+h/4; y0+h/4; y0+3*h/4; y0+3*h/4];
        rectplot([x0 x0+w-1],[y0 y0+h-1],wid);
    end
end

function rectplot(x,y,wid)
plot([x(1) x(2) x(2) x(1) x(1)],[y(1) y(1) y(2) y(2) y(1)],'c','LineWidth',wid);


function TPS  = Testing_ThinPlateSpline(ctr_ind, h, w, thrsh_sub_array)
% Thin plate spline interpolation
x = ctr_ind(:,2);
y = ctr_ind(:,1);
z = thrsh_sub_array;

%[h,w] = size(I);

nPoints=length(x);

A = zeros(nPoints);
for i=1:nPoints,
    for j=1:nPoints,
        r2 = (x(i)-x(j))^2 + (y(i)-y(j))^2;
        A(i,j) = r2*log(r2+eps);
    end
end

B = zeros(nPoints,3);
for i=1:nPoints,
    B(i,1)=1; B(i,2)=x(i); B(i,3)=y(i);
end

% find out TPS coefficients
invA = pinv(A);
b = pinv(B'*invA*B)*B'*invA*z;
a = invA*(z-B*b);

% show TPS surface
[xgrid,ygrid] = meshgrid(1:w,1:h);
R2=zeros(h,w);
TPS = b(1) + b(2)*xgrid + b(3)*ygrid;
for i=1:nPoints,
    R2 = (xgrid-x(i)).^2 + (ygrid-y(i)).^2;
    TPS(:) = TPS + a(i)*R2.*log(R2+eps);
end



% Training_GetWdthHist - Derive the histogram of stroke width from an input
% image, based on the distance transform
% 
% Usage: [vec, binloc] = Training_GetWdthHist(im)
%
% This function derives the stroke width from an input ground truth 
% images based on distance transform, and generte a histogram from the
% distance values
%
% Return: 
%              vec: the histogram of stroke width
%              binloc: the positions of bins
% Parameters: 
%              im: an grey-scale image 
%
% 
% Written by: Jichuan Shi
% Date: Jan-6-2010

function [vec, binloc] = Training_GetWdthHist(im)


%im = imread('images/Yan/gt_ocrArtificial.bmp');
dist_im = bwdist(im, 'chessboard');
sum_dist = sum(dist_im(:));
if isinf(sum_dist) || (sum_dist==0)
    vec = [];
    binloc = [];
    return;
else
    [vec, binloc] = hist_integer(dist_im(:));

    % Chop off the bin represents 0
    lgth_vec = length(vec);
    vec = vec(2:lgth_vec);
    binloc = binloc(2:lgth_vec);
    vec = vec./sum(vec);
    %figure, bar(binloc, vec, 0.1);
end

function dist_vec = Testing_HistogramDist(vec_sub, vec_gt)

% Make the vector of ground truth and segmented image have the same length
lgth_sub = length(vec_sub);
lgth_gt = length(vec_gt);

if lgth_sub-lgth_gt>0
    vec_gt = [vec_gt, zeros(1, (lgth_sub-lgth_gt))];
elseif lgth_sub-lgth_gt<0
    vec_sub = [vec_sub, zeros(1, (lgth_gt-lgth_sub))];
end

[Bhatt_dist, Hellinger_dist] = get_BhattacharyyaDis(vec_sub, vec_gt);
dist_vec = Hellinger_dist; %Bhatt_dist;

function [disBhatta disHellinger] = get_BhattacharyyaDis(seg_cur, gt_cur)

% Calculate the Bhattacharyya distance between PSD curves derived from
% segmented images and ground truth images.
%
% Input: seg_cur: the discrete PSD curve derived from segmented images
%        gt_cur:  the discrete PSD curve derived from ground truth images
% Output: disBhatta: the Bhattacharyya distance between these two PSD
%                    curves
%         disHellinger: the Hellinger distance between these two PSD
%         curves, which obeys the triangle inequality
%
% Written by: Jichuan Shi
% Date: Nov-23-2009

if sum(seg_cur)>1
% % Derive the density distribution from PSD curve
seg_dens = generate_PSD_densitycurve(seg_cur);
gt_dens = generate_PSD_densitycurve(gt_cur);
else
    seg_dens = seg_cur;
    gt_dens = gt_cur;
end

% Calculate the Bhattacharyya distance
BC = 0; %Initilize the Bhattacharyya correlation
for i = 1:length(seg_dens)
    BC = BC + (seg_dens(i)*gt_dens(i))^0.5;
end
if BC~=0
    disBhatta = -log(BC);
else
    disBhatta = inf;
end
disHellinger = (1-BC)^0.5;


% Get a histogram information from an input vector
%
% Note: the input vector contains integers only
%
% Input:  input_vec: the input vector, contains integers only
% Output: vec: the frequency counts
%         binloc: the bin location
%
% Written by: Jichuan Shi
% Date: Jan-4-2009

function [vec, binloc] = hist_integer(input_vec)

nbin = max(input_vec); % Get the maximize number of input vector as the number of bins

vec = zeros(1, nbin+1);
for i = 1:length(input_vec)
    vec(input_vec(i)+1) = vec(input_vec(i)+1) + 1;
end

binloc = zeros(1, nbin+1);
for i = 2:(nbin+1)
    binloc(i) = i-1;
end

function densCurve = generate_PSD_densitycurve(psdCurve)
% Get the probability density values from Particle Size Distribution
% Curve(PSD), because the PSD curve is an accumulative probability
% ditribution.
%
% Input: psdCurve: PSD curve which is an accumulative distribution
% Output: densCurve: the probability density distribution curve derived
% from the PSD
%
% Written by: Jichuan Shi
% Date Nov-23-2009

densCurve = zeros(size(psdCurve));
for i = 2:length(psdCurve)
    densCurve(i) = psdCurve(i) - psdCurve(i-1);    
end