-
Notifications
You must be signed in to change notification settings - Fork 6
/
hausdorff.m
132 lines (120 loc) · 3.23 KB
/
hausdorff.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
function [D,idx] = hausdorff(A,B,method)
%HAUSDORFF distance for point clouds.
%
%function [D,IDX]=hausdorff(a,b,method='euclidean')
%
% D = HAUSDORFF(A,B) computes the Hausdorff distance between
% point sets A and B. Rows of A and B correspond to observations,
% and columns correspond to variables. A and B must have same number
% of columns.
%
% D = HAUSDORFF(A,B,METHOD) lets you compute the Hausdorff distance
% with an alternate point-to-point distance. METHOD can be any
% method supported by PDIST2. METHOD defaults to 'euclidean' if not
% specified.
%
% D = HAUSDORFF(A,B,DISTFUN) lets you compute the Hausdorff distance
% with a distance function specified using a function handle @
%
% [D,IDX] = HAUSDORFF(...) also returns the indices of the farthest
% points contributing to the distance.
%
% Notes
% -----
% HAUSDORFF uses PDIST2 for computation. For gridded image data
% it is often preferred to use IMHAUSDORFF.
%
% HAUSDORFF(A,[]) = inf
% HAUSDORFF([],[]) = 0
%
% Example 1
% ---------
% Compute the Hausdorff distance for a binary segmentation.
% This is much slower than using IMHAUSDORFF.
%
% % Read in an image with an object we wish to segment.
% A = imread('hands1.jpg');
%
% % Convert the image to grayscale.
% I = rgb2gray(A);
%
% % Use active contours to segment the hand.
% mask = false(size(I));
% mask(25:end-25,25:end-25) = true;
% BW = activecontour(I, mask, 300);
%
% % Read in the ground truth against which to compare the segmentation.
% BW_groundTruth = imread('hands1-mask.png');
%
% % Extract object point coordinates
% A=regionprops(BW,'PixelList');
% B=regionprops(BW_groundTruth,'PixelList');
%
% % Compute the Hausdorff distance of this segmentation.
% [distance,idx] = hausdorff(A.PixelList,B.PixelList);
%
% % Display both masks on top of one another.
% figure
% imshowpair(BW, BW_groundTruth)
% title(['Hausdorff distance = ' num2str(distance)])
%
% % Display a line indicating the farthest points
% p=[A.PixelList(idx(1),:);B.PixelList(idx(2),:)];
% hold on;
% plot(p(:,1),p(:,2),'rx-','linewidth',2');
% hold off
%
% See also IMHAUSDORFF, PDIST2.
%
%Author: Joakim Lindblad
% Copyright (c) 2019, Joakim Lindblad
if size(A,2) ~= size(B,2)
error('A and B must have the same number of columns.');
end
if nargin < 3
method = 'euclidean';
else
if strcmp(method,'chessboard')
method = 'chebychev'; % synonymous
end
end
if isempty(A) || isempty(B)
if isempty(A) && isempty(B)
HD=0;
else
HD=inf;
end
return
end
if strcmp(method,'euclidean')
method='squaredeuclidean'; % faster
apply_root=true;
else
apply_root=false;
end
% Max of dist from A to B and B to A
if (size(A,1)*size(B,1) < 1e8)
D = pdist2(A,B,method);
[D1,idxA1] = min(D,[],1);
[D2,idxB1] = min(D,[],2);
clear D;
else
% Less memory hungry version
[D1,idxA1] = pdist2(A,B,method,'Smallest',1);
[D2,idxB1] = pdist2(B,A,method,'Smallest',1);
end
[D1,idxB2]=max(D1);
[D2,idxA2]=max(D2);
if (D1>D2)
D=D1;
idx=[idxA1(idxB2),idxB2];
else
D=D2;
idx=[idxA2,idxB1(idxA2)];
end
if apply_root
D = sqrt(D);
end
if (nargout < 2)
clear idx;
end