source: prextra/hausdm.m @ 5

%HAUSDM Hausdorff distance between datasets of image blobs
%
%       [DH,DM] = hausdm(A,B,fid)
%
% Computes a Hausdorff distance matrix between the sets of binary images 
% A and B, or datasets containing them as features. 
% Progress is reported in fid (fid = 1: on the sreeen).
% If A and B are image sets and 
% size(A) is [may,max,na]
% size(B) is [mby,mbx,nb]
% then DH is the Hausdorff distance matrix, size(DH) = [na,nb], between 
% these sets. DM is the Modified Hausdorff distance matrix.
% Preferably na <= nb (faster).
%
% See M.-P. Dubuisson and A.K. Jain, Modified Hausdorff distance for object
%     matching, Proceedings 12th IAPR International Conference on Pattern
%     Recognition (Jerusalem, October 9-13, 1994), vol. 1, IEEE, Piscataway,
%     NJ, USA,94CH3440-5, 1994, 566-568.
%
% $Id: hausdm.m,v 1.1 2005/04/25 05:54:36 duin Exp $

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function [dh,dm] = hausdm(A,B,fid)
if nargin < 3, fid = 0; end

if isdataset(A) & isdataset(B)
        [dh,dm] = hausdm(data2im(A),data2im(B),fid);
        dh = setdata(A,dh);
        dm = setdata(A,dm);
        return
end

[ma1,ma2,na] = size(A);
[mb1,mb2,nb] = size(B);
dh = zeros(na,nb);
dm = zeros(na,nb);
for i=1:na
        a = A(:,:,i); 
        J = find(any(a));
        J = [min(J):max(J)];
        K = find(any(a'));
        K = [min(K):max(K)];
        a = double(a(K,J));
        if length(a(:)) > 0
                a = bord(a,0);
        end
        ca = contourc(a,[0.5,0.5]);
        J = find(ca(1,:) == 0.5);
        ca(:,[J J+1]) =[];
        ca = ca - repmat([1.5;1.5],1,size(ca,2));
        ca = ca/max(ca(:));
        ca = ca - repmat(max(ca,[],2)/2,1,size(ca,2));
        for j = 1:nb
                b = B(:,:,j);
                J = find(any(b));
                J = [min(J):max(J)];
                K = find(any(b'));
                K = [min(K):max(K)];
                b = double(b(K,J));
                if length(b(:)) > 0
                        b = bord(b,0);
                end
                cb = contourc(b,[0.5,0.5]);
                J = find(cb(1,:) == 0.5);
                cb(:,[J J+1]) =[];
                cb = cb - repmat([1.5;1.5],1,size(cb,2));
                cb = cb/max(cb(:));
                cb = cb - repmat(max(cb,[],2)/2,1,size(cb,2));
                dab = sqrt(distm(ca',cb'));
                dh(i,j) = max(max(min(dab)),max(min(dab')));
                dm(i,j) = max(mean(min(dab)),mean(min(dab')));
%               if fid, disp([i,j,dh(i,j),dm(i,j)]); end
                fprintf(fid,'%5d %5d %10.3f %8.3f \n',i,j,dh(i,j),dm(i,j));
        end
end
        
% C = bord(A,n,m)
% Puts a border of width m (default m=1) around image A
% and gives it value n. If n = NaN: mirror image values.
function C = bord(A,n,m);
%ipcontr(0);
if nargin == 2; m=1; end
[x,y] = size(A);
if m > min(x,y)
        mm = min(x,y);
        C = bord(A,n,mm);
        C = bord(C,n,m-mm);
        return
end
if isnan(n)
   C = [A(:,m:-1:1),A,A(:,y:-1:y-m+1)];
   C = [C(m:-1:1,:);C;C(x:-1:x-m+1,:)];
else
   bx = ones(x,m)*n;
   by = ones(m,y+2*m)*n;
   C = [by;[bx,A,bx];by];
end
return