%Written by: Sukanta Basu (Last updated: Jan 14, 2007)

%**************************************************************************
%Computes Partial Hausdorff distances between 2 images A & B (binary)
%In our case: A=Observed Rainfall, B=Forecasted Rainfall; Also,
%computes the denominator (a modified version of UIQI - see below) in
%FQI.
%**************************************************************************

function [dH,dQ1,dQ1x] = Hausdorff(A,B,prc)

%**************************************************************************
% dH - Partial Hausdorff Distance (Numerator) - Distance measure/error.
% dQ1 - UIQI (Denominator)                    - Amplitude measure/error
%                                             - Uses Non-Zero pixels.
% dQ2 - UIQI (Denominator)                    - Amplitude measure/error.
%                                             - Uses All pixels.
%**************************************************************************

[ai,aj,av]=find(A); Na=length(av);
[bi,bj,bv]=find(B); Nb=length(bv);

%Forward Distance**********************************************************
if Na*Nb < 1e6
    df = abs(ones(Nb,1)*ai' - bi*ones(1,Na)) + abs(ones(Nb,1)*aj' - bj*ones(1,Na));
    df = min(df,[],1);
else
    df = inf*ones(Na,1); d = zeros(Na,1);
    for k1 = 1:Na
        for k2 = 1:Nb
            d(k1) = abs(ai(k1)-bi(k2)) + abs(aj(k1)-bj(k2));
            if d(k1) < df(k1)
                df(k1) = d(k1);
            end
        end
    end
end
df = sort(df);
Indx = round(prc*length(df));                   
df = df(Indx);
%**************************************************************************

%Backward Distance*********************************************************
if Na*Nb < 1e6
    db = abs(ones(Na,1)*bi' - ai*ones(1,Nb)) + abs(ones(Na,1)*bj' - aj*ones(1,Nb));
    db = min(db,[],1); 
else
    db = inf*ones(Nb,1); d = zeros(Nb,1);
    for k1 = 1:Nb
        for k2 = 1:Na
            d(k1) = abs(bi(k1)-ai(k2)) + abs(bj(k1)-aj(k2));
            if d(k1) < db(k1)
                db(k1) = d(k1);
            end
        end
    end    
end
db   = sort(db);	
Indx = round(prc*length(db));						
db   = db(Indx); 
dH   = max(df,db);                                  %Hausdorff Distance
%**************************************************************************

%**************************************************************************
% Amplitude Error
% Computes modified Universal Image Quality Index (UIQI), which appears as the
% denominator in FQI
% see Wang, Z., and A. C. Bovik, A universal image quality index,
% IEEE Signal Processing Letters,  9, 81-84, 2002. 
%**************************************************************************
% Similar to JGR Paper: Only uses Non-Zero Pixels
mA  = mean(av);                                             %Non-zero pixel mean
mB  = mean(bv);
sigA= sqrt( (1/(Na-1)) * (sum((av-mA).^2)));				%Non-zero pixel stdev
sigB= sqrt( (1/(Nb-1)) * (sum((bv-mB).^2)));
F2  = 2*mA*mB/(mA^2+mB^2);
F3  = 2*sigA*sigB/(sigA^2+sigB^2);
dQ1 = (F2*F3);

% Use 'All' Pixels
mAx  = mean(mean(A));
mBx  = mean(mean(B));
sigAx= sqrt( (1/(numel(A)-1))*sum(sum((A-mAx).^2)));
sigBx= sqrt( (1/(numel(B)-1))*sum(sum((B-mBx).^2)));
F2x  = 2*mAx*mBx/(mAx^2+mBx^2);
F3x  = 2*sigAx*sigBx/(sigAx^2+sigBx^2);
dQ1x = (F2x*F3x);
%**********************************************************************************