豪斯多夫距离-- Hausdorff distance of convex polygons

蒙特利尔的麦吉尔大学的计算几何课程资料：

原文链接：http://cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/98/normand/main.html

 

1.  Introduction

When talking about distances, we usually mean the shortest :   for instance, if a point X is said to be at distance D of a polygon P, we generally assume that D is the distance from X to the nearest point of P. The same logic applies for polygons :  if two polygons A and B are at some distance from each other, we commonly understand that distance as the shortest one between any point of A and any point of B.  Formally, this is called a minimin function, because the distance D between A and B is given by :

(eq. 1)

This equation reads like a computer program : « for every point a of A, find its smallest distance to any point b of B ;  finally, keep the smallest distance found among all points a ».

That definition of distance between polygons can become quite unsatisfactory for some applications ;  let's see for example fig. 1.  We could say the triangles are close to each other considering their shortest distance, shown by their red vertices. However, we would naturally expect that a small distance between these polygons means that no point of one polygon is far from the other polygon.   In this sense, the two polygons shown in fig. 1 are not so close, as their furthest points, shown in blue, could actually be very far away from the other polygon. Clearly, the shortest distance is totally independent of each polygonal shape.

Figure 1 :

The shortest distance doesn't consider the whole shape.

Another example is given by fig. 2, where we have the same two triangles at the same shortest distance than in fig. 1, but in different position.  It's quite obvious that the shortest distance concept carries very low informative content, as the distance value did not change from the previous case, while something did change with the objects.

Figure 2 :

The shortest distance doesn't account for the position of the objects.

As we'll see in the next section, in spite of its apparent complexity, the Hausdorff distance does capture these subtleties, ignored by the shortest distance.

2.  What is Hausdorff distance ?

Named after Felix Hausdorff (1868-1942), Hausdorff distance is the «  maximum distance of a set to the nearest point in the other set » [Rote91]. More formally, Hausdorff distance from set A to set B is a maximin function, defined as

(eq. 2)

where a and b are points of sets A and B respectively, and d(a, b) is any metric between these points ; for simplicity, we'll take d(a, b) as the Euclidian distance between a and b. If for instance A and B are two sets of points, a brute force algorithm would be :

Brute force algorithm : 1.  h = 0  2.  for every point ai of A,       2.1  shortest = Inf ;       2.2  for every point bj of B                     dij = d (ai , bj )                     if dij < shortest then                               shortest = dij       2.3  if shortest > h then                      h = shortest

Figure 3 :   Hausdorff distance on point sets.

This is illustrated in fig. 3 : just click on the arrow to see the basic steps of this computation. This algorithm obviously runs in O(n m) time, with n and m the number of points in each set. 

It should be noted that Hausdorff distance is oriented (we could say asymmetric as well), which means that most of times h(A, B) is not equal to h(B, A). This general condition also holds for the example of fig. 3, as h(A, B) = d(a1, b1), while h(B, A) = d(b2, a1). This asymmetry is a property of maximin functions, while minimin functions are symmetric.

A more general definition of Hausdorff distance would be :

H (A, B) = max { h (A, B), h (B, A) }

(eq. 3)

which defines the Hausdorff distance between A and B, while eq. 2 applied to Hausdorff distance from A to B (also called directed Hausdorff distance). The two distances h(A, B) and h(B, A) are sometimes termed as forward and backward Hausdorff distances of A to B. Although the terminology is not stable yet among authors, eq. 3 is usually meant when talking about Hausdorff distance. Unless otherwise mentionned, from now on we will also refer to eq. 3 when saying "Hausdorff distance".

If sets A and B are made of lines or polygons instead of single points, then H(A, B) applies to all defining points of these lines or polygons, and not only to their vertices. The brute force algorithm could no longer be used for computing Hausdorff distance between such sets, as they involve an infinite number of points.

 

 

So, what about the polygons of fig. 1 ? Remember, some of their points were close, but not all of them. Hausdorff distance gives an interesting measure of their mutual proximity, by indicating the maximal distance between any point of one polygon to the other polygon. Better than the shortest distance, which applied only to one point of each polygon, irrespective of all other points of the polygons.

 

Figure 4 :

Hausdorff distance shown around extremum of each triangles of fig. 1. Each circle has a radius of H( P1 , P2).

The other concern was the insensitivity of the shortest distance to the position of the polygons. We saw that this distance doesn't consider at all the disposition of the polygons.  Here again, Hausdorff distance has the advantage of being sensitive to position, as shown in fig.5. 

 

Figure 5 :

Hausdorff distance for the triangles of fig. 4 at the same shortest distance, but in different position.

3.  Computing Hausdorff distance between convex polygons

3.1  Assumptions

Throughout the rest of our discussion, we assume the following facts about polygons A and B :

• Polygons A and B are simple convex polygons ;
   
• Polygons A and B are disjoint from each other, that is :
  
    - they don't intersect together ; 
    - no polygon contains the other.
   

 

3.2  Lemmas

The algorithm explained in the next section is based on three geometric observations, presented here. In order to simplify the text, we assume two points a and b that belong respectively to polygons A and B, such that :
 

d (a, b) = h (A, B)
 

In simple words, a is the furthest point of polygon A relative to polygon B, while b is the closest point of polygon B relative to polygon A.

 

Lemma 1a :

The perpendicular to ab at a is a supporting line of A, and A is on the same side as B relative to that line. Proof of lemma 1a

Lemma 1b :

The perpendicular to ab at b is a supporting line of B, and a and B are on different sides relative to that line. Proof of lemma 1b

Lemma 2 :

There is a vertex x of A such that the distance from x to B is equal to h (A, B). Proof of lemma 2

Lemma 3 :

Let bi be the closest point of B from a vertex a i  of A.  If µ is the moving direction (clockwise or counterclockwise) from bi to bi+1 then, for a complete cycle through all vertices of A, µ changes no more than twice. Proof of lemma 3

3.3  Algorithm

The algorithm presented here was proposed by [Atallah83].  Its basic strategy is to compute successively h(A,B) and h(B, A) ;  because of lemma 2, there is no need to query every point of the starting polygon, but only its vertices. 
 

An important fact used by this algorithm is that a closest point can only be a vertex of the target polygon, or the foot z of a line perpendicular to one of its edges.  This fact suggests a function to check for the existence of a possible closest point. Given a source point a and a target edge defined by a point b1 and a vertex b2 :

Function z = CheckForClosePoint (a, b1 , b2 ) :

Compute the position z where the line that passes through b1 and b2 crosses its perpendicular through a  ;
if z is between b1 b2 then return z ;
else compute at b2 a line P perpendicular to the line ab2 ;
        if P is a supporting line of B then return b2 ;
        else return NULL. 
 

That function obviously uses lemma 1b to decide whether or not the closest point of B might be located on the target edge, that should be close to a.  It also supposes that the source point a and b2 are not located on different sides of the perpendicular to [b1b2 ] at b1, accordingly to lemma 3. 

Now we are ready for the main algorithm ; the vertices of both polygons are presumed to be enumerated counterclockwise :

 

Algorithm for computing h(A, B) : 1.  From a1, find the closest point b1 and compute d1 = d ( a1, b1 ) 2.  h(A, B) = d1 3.  for each vertex ai of A,       3.1  if ai+1 is to the left of aibi                      find bi+1 , scanning B counterclockwise with CheckForClosePoint from bi               if ai+1 is to the right of aibi                      find bi+1 , scanning B clockwise with CheckForClosePoint from bi               if ai+1 is anywhere on aibi                       bi+1 = bi       3.2  Compute di+1 = d (ai+1 , bi+1 )       3.3  h (A, B) = max { h (A, B), di+1 }

 

3.4  Complexity

If polygons A and B respectively have n and m vertices, then :

 

• Step 1 can clearly be done in O(m) time ;
• Step 2 takes constant time O(1) ;
• Step 3 will be executed (n-1) times, that is O(n) ;
• Step 3.1 will not be executed in total more than O(2m). This is a consequence of lemma 3, which guarantees that polygon B can not be scanned more than twice ;
• Steps 3.2 and 3.3 are done in constant time O(1) ;

So the algorithm for computing h(A, B) takes :
 

O(m) + O(n) + O(2m) = O(n+m)

To find H(A, B), the algorithm needs to executed twice ;  the total complexity for computing Hausdorff distance then stays linear to O(n+m).

 

3.5  Interactive applet

This applet illustrates the algorithm for computing h(A,B). You only need to draw two polygons, and then press the "step" or "run" button. Left click to define a new vertex, and close the polygon by clicking near the first vertex. Polygon A is the first one you draw, in green, while polygon B appears next, in red. 

The algorithm was slightly modified to make it more appealing visually. Even if this algorithm is intended for two polygons totally separated from each other, it also works when B is inside A. However, it won't work if A is inside of B, or when A and B are partially intersecting. You're allowed anyway to try these cases to see what happens ! 

When defining your polygons, you will see a yellow area that indicates where you can add the next vertex, so the polygon keeps convex. The applet won't let you define a non-convex polygon. 
 

Please notice that the first time you draw the second half of a polygon, you will have to wait a few seconds until the Jama package loads. 

 

 

 

 

 

 

 

1. import java.awt.*;

2. import java.applet.Applet;

3. import java.util.Vector;

4. import java.lang.Math;

5. import Jama.*;

6.  

7. public class Hausdorff extends Applet

8. {

9. PolyArea area;

10. Panel control;

11. Button runStop;

12. boolean running;

13. TextField comment;

14.  

15. public void init()

16. {

17. setLayout (new BorderLayout());

18. comment = new TextField ("", 60);

19. comment.setEditable (false);

20. area = new PolyArea (comment);

21. add ("Center", area);

22.  

23. control = new Panel();

24. control.add (new Button ("step"));

25. runStop = new Button ("run");

26. control.add (runStop);

27. control.add (new Button ("reset"));

28. control.add (comment);

29. add("South", control);

30.  

31. running = false;

32. }

33.  

34. public boolean action (Event evt, Object arg)

35. {

36. if ("step".equals(arg))

37. area.stepAlgo();

38.  

39. if ("run".equals(arg))

40. {

41. runStop.setLabel ("stop");

42. startAnim();

43. running = true;

44. }

45.  

46. if ("stop".equals(arg))

47. {

48. runStop.setLabel ("run");

49. stopAnim();

50. running = false;

51. }

52.  

53. if ("reset".equals(arg))

54. {

55. remove (area);

56. stopAnim();

57. area = new PolyArea (comment);

58. add ("Center", area);

59. runStop.setLabel ("run");

60. validate();

61.  

62. stopAnim();

63. running = false;

64. }

65.  

66. return true;

67. }

68.  

69. public void start()

70. {

71. if (running)

72. startAnim();

73. }

74.  

75. public void stop()

76. {

77. stopAnim();

78. }

79.  

80. public void startAnim()

81. {

82. if (area.animator == null)

83. {

84. area.animator = new Thread (area);

85. }

86. area.animator.start();

87. }

88.  

89. public void stopAnim()

90. {

91. area.animator = null;

92. }

93.  

94. public void paint(Graphics g) {

95. Dimension d = getSize();

96. g.setColor (Color.black);

97. g.drawRect(0,0, d.width - 1, d.height - 1);

98. }

99.  

100. /* This is used to leave room to the black box painted in

101. * the paint() method. If we don't do that, it is overwritten.

102. */

103. public Insets getInsets() {

104. return new Insets(3,3,3,3);

105. }

106. }

107.  
108.  

109. //=============================================================================

110.  
111.  

112. /*

113. * The PolyArea class defines an area that will hold our two polygons.

114. * It will first create them by catching mouse clicks events and adding

115. * the points to the polygons, and it will then run the algorithm on the

116. * polygons.

117. */

118. class PolyArea extends Canvas implements Runnable

119. {

120. Dimension offDimension;

121. Image offImage;

122. Graphics offGraphics;

123.  

124. TextField comment;

125.  

126. public Thread animator;

127.  

128. FConvexPoly poly1, poly2;

129.  

130. int nextStep;

131. int currentV1;

132. FPoint bestV1, bestV2, currentV2;

133. double bestLength;

134. int currentRegBase;

135. boolean band;

136. Polygon currentRegion;

137. boolean trigo;

138.  

139. public PolyArea (TextField comment)

140. {

141. animator = null;

142. this.comment = comment;

143.  

144. setBackground (Color.white);

145.  

146. poly1 = new FConvexPoly (Color.green);

147. poly2 = new FConvexPoly (Color.red);

148.  

149. nextStep = -1;

150. currentV1 = currentRegBase = -1;

151. bestV1 = bestV2 = currentV2 = null;

152. band = false;

153. currentRegion = null;

154. bestLength = 0;

155. trigo = true;

156.  

157. comment.setText ("Please enter the first polygon");

158. comment.setBackground (new Color (220, 255, 230));

159. comment.setForeground (Color.black);

160. }

161.  

162. public boolean handleEvent (Event e)

163. {

164. switch (e.id)

165. {

166. case Event.MOUSE_DOWN:

167. if (nextStep == -1)

168. {

169. if (! poly1.isClosed())

170. {

171. poly1.addPoint (new FPoint (e.x, e.y));

172. if (poly1.isClosed())

173. {

174. comment.setText ("Please enter the second polygon");

175. comment.setBackground (new Color (255, 220, 220));

176. }

177. }

178. else

179. {

180. poly2.addPoint (new FPoint (e.x, e.y));

181. if (poly2.isClosed())

182. {

183. comment.setText ("Press step or run to see the algorithm");

184. comment.setBackground (new Color (255, 255, 220));

185. nextStep = 0;

186. }

187. }

188. }

189.  

190. repaint();

191. return true;

192.  

193. default:

194. return false;

195. }

196. }

197.  

198. /* This method performs a step in the algorithm. It is called either

199. * by the applet when the "step" button is clicked, or by the animator

200. * thread if this one is running. */

201. public void stepAlgo ()

202. {

203. Point p;

204.  

205. switch (nextStep)

206. {

207. case 0:

208. comment.setText ("Searching for the first vertex");

209. comment.setBackground (new Color (255, 235, 200));

210.  

211. currentV1 = 0;

212. currentRegBase = 0;

213. band = false;

214. makeRegion();

215.  

216. p = poly1.getPoint (currentV1);

217. if (currentRegion.inside (p.x, p.y))

218. nextStep = 2;

219. else

220. nextStep = 1;

221.  

222. repaint();

223. break;

224.  

225. case 1:

226. if (! band)

227. band = true;

228. else

229. {

230. currentRegBase++;

231. band = false;

232. }

233. makeRegion();

234. p = poly1.getPoint (currentV1);

235. if (currentRegion.inside (p.x, p.y))

236. nextStep = 2;

237. else

238. nextStep = 1;

239.  

240. repaint();

241. break;

242.  

243. case 2:

244. comment.setText ("Computing length");

245. comment.setBackground (new Color (255, 220, 255));

246.  

247. makeV2();

248. bestV1 = poly1.getFPoint (currentV1);

249. bestV2 = currentV2;

250. bestLength = new FVector (bestV1, bestV2).length();

251.  

252. nextStep = 3;

253. repaint();

254. break;

255.  

256. case 3:

257. comment.setText ("Searching for the next vertex");

258. comment.setBackground (new Color (255, 235, 200));

259.  

260. if (new FVector (currentV2, poly1.getFPoint (currentV1)).isTrigo

261. (new FVector (currentV2, poly1.getFPoint (currentV1 + 1))))

262. trigo = true;

263. else

264. trigo = false;

265.  

266. currentV1++;

267. currentV2 = null;

268. if (poly1.getFPoint (currentV1) == poly1.getFPoint (0))

269. nextStep = 7;

270. else

271. {

272. p = poly1.getPoint (currentV1);

273. if (currentRegion.inside (p.x, p.y))

274. nextStep = 5;

275. else

276. nextStep = 4;

277. }

278.  

279. repaint();

280. break;

281.  

282. case 4:

283. if ((trigo || !poly2.isTrigo()) && !(trigo && !poly2.isTrigo()))

284. if (! band)

285. band = true;

286. else

287. {

288. currentRegBase++;

289. band = false;

290. }

291. else

292. if (! band)

293. {

294. currentRegBase--;

295. band = true;

296. }

297. else

298. band = false;

299.  

300. makeRegion();

301. p = poly1.getPoint (currentV1);

302. if (currentRegion.inside (p.x, p.y))

303. nextStep = 5;

304. else

305. nextStep = 4;

306.  

307. repaint();

308. break;

309.  

310. case 5:

311. comment.setText ("Comparing this length with the previous best");

312. comment.setBackground (new Color (255, 220, 255));

313.  

314. makeV2();

315.  

316. if (new FVector (poly1.getFPoint (currentV1), currentV2).length() > bestLength)

317. nextStep = 6;

318. else

319. nextStep = 3;

320.  

321. repaint();

322. break;

323.  

324. case 6:

325. comment.setText ("This is the new best length");

326. comment.setBackground (new Color (220, 255, 220));

327.  

328. bestV1 = poly1.getFPoint (currentV1);

329. bestV2 = currentV2;

330. bestLength = new FVector (bestV1, bestV2).length();

331.  

332. nextStep = 3;

333. repaint();

334. break;

335.  

336. case 7:

337. comment.setText ("Hausdorff distance from poly 1 to poly 2");

338. comment.setBackground (new Color (220, 220, 255));

339.  

340. currentV1 = -1;

341. currentV2 = null;

342. currentRegion = null;

343. animator = null;

344.  

345. repaint();

346. break;

347.  

348. default:

349. break;

350. }

351. }

352.  

353. /* The paint method draws the current state of the algorithm in the

354. * given Graphics, including the two polygons, the yellow region and

355. * the important points. */

356. public void paint (Graphics g)

357. {

358. poly1.fill (g);

359. poly2.fill (g);

360.  

361. if (currentRegion != null)

362. {

363. g.setColor (Color.yellow);

364. g.fillPolygon (currentRegion);

365. }

366.  

367. poly1.draw (g);

368. poly2.draw (g);

369.  

370. if (currentV1 >= 0)

371. {

372. Point p = poly1.getPoint (currentV1);

373. g.setColor (Color.blue);

374. g.drawOval (p.x-5, p.y-5, 10, 10);

375. }

376.  

377. if (currentV2 != null)

378. {

379. Point p1 = poly1.getPoint (currentV1);

380. Point p2 = currentV2.getPoint();

381. g.setColor (Color.blue);

382. g.drawOval (p2.x-5, p2.y-5, 10, 10);

383. g.setColor (Color.blue);

384. g.drawLine (p1.x, p1.y, p2.x, p2.y);

385. }

386.  

387. if (bestV1 != null)

388. {

389. Point p1 = bestV1.getPoint();

390. Point p2 = bestV2.getPoint();

391. g.setColor (Color.magenta);

392. g.drawLine (p1.x, p1.y, p2.x, p2.y);

393. }

394. }

395.  

396. /* When called by the AWT, the update method clears its offscreen

397. * buffer, calls paint() to draw in it, and then copies the buffer to the

398. * screen. */

399. public void update (Graphics g)

400. {

401. Dimension d = size();

402.  

403. if ((offGraphics == null) ||

404. (d.width != offDimension.width) ||

405. (d.height != offDimension.height))

406. {

407. offDimension = d;

408. offImage = createImage(d.width, d.height);

409. offGraphics = offImage.getGraphics();

410. }

411.  

412. offGraphics.setColor (getBackground());

413. offGraphics.fillRect (0, 0, d.width, d.height);

414. paint (offGraphics);

415.  

416. g.drawImage (offImage, 0, 0, this);

417. }

418.  

419. /* This method is called in the algorithm to make the new

420. * yellow polygon that corresponds to the sweeping yellow

421. * region. */

422. void makeRegion()

423. {

424. Polygon region;

425. FPoint p1, p2;

426. FVector v1, v2, edge;

427.  

428. p1 = poly2.getFPoint (currentRegBase);

429. p2 = poly2.getFPoint (currentRegBase + 1);

430. edge = new FVector (p1, p2);

431.  

432. if (poly2.isTrigo())

433. v1 = edge.normal().mult(-1);

434. else

435. v1 = edge.normal();

436.  

437. if (band)

438. {

439. currentRegion = FConvexPoly.infiniteRegion (p1, v1, p2, v1);

440. }

441. else

442. {

443. p2 = poly2.getFPoint (currentRegBase - 1);

444. edge = new FVector (p2, p1);

445. if (poly2.isTrigo())

446. v2 = edge.normal().mult(-1);

447. else

448. v2 = edge.normal();

449.  

450. currentRegion = FConvexPoly.infiniteRegion (p1, v1, p1, v2);

451. }

452. }

453.  

454. /* This method creates the new end of the current smallest distance

455. * that's on the red polygon. If the current vertex of the green polygon

456. * is in a yellow sector, it corresponds to the current vertex of the red

457. * polygon. If it is in a band, it is the foot of the perpendicular to the

458. * current edge of the red polygon that goes through the vertex. */

459. void makeV2()

460. {

461. if (! band)

462. currentV2 = poly2.getFPoint (currentRegBase);

463. else

464. {

465. FPoint p1, p2;

466. FVector vBase;

467. FLine baseLine, perpendicular;

468.  

469. p1 = poly2.getFPoint (currentRegBase);

470. p2 = poly2.getFPoint (currentRegBase + 1);

471. vBase = new FVector (p1, p2);

472. baseLine = new FLine (p1, vBase);

473.  

474. p2 = poly1.getFPoint (currentV1);

475. perpendicular = new FLine (p2, vBase.normal());

476.  

477. currentV2 = baseLine.intersect (perpendicular);

478. }

479. }

480.  

481. /* The run method is called by the virtual machine to perform the

482. * automated animation of the algorithm. It calls the stepAlgo() method

483. * three times per second to animate the algorithme. */

484. public void run()

485. {

486. Thread.currentThread().setPriority(Thread.MIN_PRIORITY);

487.  

488. long startTime = System.currentTimeMillis();

489.  

490. while (Thread.currentThread() == animator)

491. {

492. stepAlgo();

493.  

494. try

495. {

496. startTime += 300;

497. Thread.sleep (Math.max (0, startTime-System.currentTimeMillis()));

498. }

499. catch (InterruptedException e) { break; }

500. }

501. }

502. }

503.  
504.  

505. //=============================================================================

506.  
507.  

508. /*

509. * The FConvexPoly class represents a convex polygon.

510. */

511. class FConvexPoly

512. {

513. public Vector v;

514. boolean closed;

515. boolean trigo;

516. Color polyColor;

517.  

518. public FConvexPoly (Color c)

519. {

520. v = new Vector();

521. closed = false;

522. trigo = false;

523. polyColor = c;

524. }

525.  

526. public boolean isClosed ()

527. {

528. return closed;

529. }

530.  

531. public boolean isTrigo ()

532. {

533. return trigo;

534. }

535.  

536. public Point getPoint (int i)

537. {

538. return ((FPoint) v.elementAt (MesMath.mod (i, v.size()))).getPoint();

539. }

540.  

541. public FPoint getFPoint (int i)

542. {

543. return (FPoint) v.elementAt (MesMath.mod (i, v.size()));

544. }

545.  

546. public void addPoint (FPoint p)

547. {

548. if (! closed)

549. {

550. if (v.size() < 3)

551. {

552. v.addElement (p);

553. if (v.size() == 3)

554. trigo = firstVector().isTrigo (lastVector());

555. }

556. else

557. {

558. if (p.equals ((FPoint) v.firstElement()))

559. closed = true;

560. else

561. {

562. if (isNextOK (p))

563. v.addElement (p);

564. }

565. }

566. }

567. }

568.  

569. public boolean isNextOK (FPoint p)

570. {

571. FVector next;

572.  

573. next = new FVector ((FPoint) v.firstElement(), p);

574. if (firstVector().isTrigo (next) != trigo)

575. return false;

576.  

577. next = new FVector ((FPoint) v.lastElement(), p);

578. if (lastVector().isTrigo (next) != trigo)

579. return false;

580.  

581. next = new FVector ((FPoint) v.firstElement(), p);

582. if (boundVector().isTrigo (next) == trigo)

583. return false;

584. else

585. return true;

586. }

587.  

588. public FVector firstVector()

589. {

590. return new FVector ((FPoint) v.elementAt (0),

591. (FPoint) v.elementAt (1));

592. }

593.  

594. public FVector lastVector()

595. {

596. return new FVector ((FPoint) v.elementAt (v.size()-2),

597. (FPoint) v.elementAt (v.size()-1));

598. }

599.  

600. public FVector boundVector()

601. {

602. return new FVector ((FPoint) v.elementAt (v.size()-1),

603. (FPoint) v.elementAt (0));

604. }

605.  

606. /* Draw the filled polygon on the given graphics. */

607. public void fill (Graphics g)

608. {

609. Polygon p = new Polygon();

610. Point pt, pt2;

611. FPoint fpt;

612.  

613. for (int a=0; a<v.size(); a++)

614. {

615. pt = ((FPoint) v.elementAt (a)).getPoint();

616. p.addPoint (pt.x, pt.y);

617. }

618.  

619. g.setColor (polyColor);

620. g.fillPolygon (p);

621.  

622. if (!closed && v.size() >= 3)

623. {

624. Polygon nextReg = new Polygon();

625.  

626. if (lastVector().isTrigo (firstVector()) == trigo)

627. {

628. pt = ((FPoint) v.firstElement()).getPoint();

629. pt2 = ((FPoint) v.lastElement()).getPoint();

630. nextReg.addPoint (pt.x, pt.y);

631. nextReg.addPoint (pt2.x, pt2.y);

632.  

633. FLine l1 = new FLine ((FPoint) v.firstElement(),

634. firstVector());

635. FLine l2 = new FLine ((FPoint) v.lastElement(),

636. lastVector());

637.  

638. fpt = l1.intersect (l2);

639.  

640. pt = fpt.getPoint();

641. nextReg.addPoint (pt.x, pt.y);

642. }

643. else

644. {

645. FVector firstInv = firstVector().mult (-1);

646. nextReg = infiniteRegion ((FPoint) v.firstElement(),

647. firstInv,

648. (FPoint) v.lastElement(),

649. lastVector());

650. }

651.  

652. g.setColor (Color.yellow);

653. g.fillPolygon (nextReg);

654. }

655. }

656.  

657. /* Draw the outline of the polygon on the given graphics. */

658. public void draw (Graphics g)

659. {

660. Point pt, pt2;

661.  

662. g.setColor (Color.blue);

663. if (v.size() > 0)

664. {

665. pt = pt2 = ((FPoint) v.firstElement()).getPoint();

666.  

667. for (int a=0; a<v.size()-1; a++)

668. {

669. pt2 = ((FPoint) v.elementAt (a+1)).getPoint();

670. g.drawLine (pt.x, pt.y, pt2.x, pt2.y);

671. pt = pt2;

672. }

673. pt = ((FPoint) v.firstElement()).getPoint();

674.  

675. if (closed)

676. g.drawLine (pt.x, pt.y, pt2.x, pt2.y);

677. else

678. {

679. g.setColor (Color.red);

680. g.drawOval (pt.x-5, pt.y-5, 10, 10);

681. }

682. }

683. }

684.  

685. /* This method returns a Polygon that will look infinite when

686. * drawn by the Graphics methods. This should not be in here, it should

687. * not be a static method, but it worked so I left it this way...*/

688. public static Polygon infiniteRegion (FPoint p1, FVector v1, FPoint p2, FVector v2)

689. {

690. Polygon p = new Polygon();

691. Point pt;

692. double length;

693. FVector middle;

694.  

695. pt = p1.getPoint();

696. p.addPoint (pt.x, pt.y);

697.  

698. if (! p1.equals (p2))

699. {

700. pt = p2.getPoint();

701. p.addPoint (pt.x, pt.y);

702. }

703.  

704. length = v2.length();

705. pt = new Point ((int) (p2.x + (v2.x / length * 2000)),

706. (int) (p2.y + (v2.y / length * 2000)));

707. p.addPoint (pt.x, pt.y);

708.  

709. middle = (v1.mult (1/v1.length())).add (v2.mult (1/v2.length()));

710. length = middle.length();

711. pt = new Point ((int) (middle.x + (middle.x / length * 2000)),

712. (int) (middle.y + (middle.y / length * 2000)));

713. p.addPoint (pt.x, pt.y);

714.  

715. length = v1.length();

716. pt = new Point ((int) (p1.x + (v1.x / length * 2000)),

717. (int) (p1.y + (v1.y / length * 2000)));

718. p.addPoint (pt.x, pt.y);

719.  

720. return p;

721. }

722. }

723.  
724.  

725. //=============================================================================

726.  
727.  

728. /*

729. * The FPoint class represents a point in the plane.

730. */

731. class FPoint

732. {

733. static double EPSILON = 10;

734.  

735. double x;

736. double y;

737.  

738. public FPoint (double x, double y)

739. {

740. this.x = x;

741. this.y = y;

742. }

743.  

744. public FPoint (int x, int y)

745. {

746. this.x = (double) x;

747. this.y = (double) y;

748. }

749.  

750. public FPoint (Point p)

751. {

752. this.x = (double) p.x;

753. this.y = (double) p.y;

754. }

755.  

756. public boolean equals (FPoint p)

757. {

758. if (new FVector (this, p).length() < EPSILON)

759. return true;

760. else

761. return false;

762. }

763.  

764. public Point getPoint()

765. {

766. return new Point ((int) x, (int) y);

767. }

768. }

769.  
770.  

771. //=============================================================================

772.  
773.  

774. /*

775. * The FVector class provides a basic 2-dimensional vector type,

776. * along with some uselful methods.

777. */

778. class FVector

779. {

780. double x;

781. double y;

782.  

783. public FVector (double x, double y)

784. {

785. this.x = x;

786. this.y = y;

787. }

788.  

789. public FVector (FPoint a, FPoint b)

790. {

791. x = b.x - a.x;

792. y = b.y - a.y;

793. }

794.  

795. /* Returns the dot product of the vector with v. */

796. public double dot (FVector v)

797. {

798. return x * v.x + y * v.y;

799. }

800.  

801. /* Returns the vector normal to this vector. */

802. public FVector normal()

803. {

804. return new FVector (-y, x);

805. }

806.  

807. /* Return true if v makes an angle with the vector in the

808. * trigonometric direction (a left turn). */

809. public boolean isTrigo (FVector v)

810. {

811. if (dot (v.normal()) <= 0)

812. return true;

813. else

814. return false;

815. }

816.  

817. public double length()

818. {

819. return Math.sqrt (x*x + y*y);

820. }

821.  

822. public FVector add (FVector v)

823. {

824. return new FVector (x + v.x, y + v.y);

825. }

826.  

827. public FVector mult (double a)

828. {

829. return new FVector (x*a, y*a);

830. }

831. }

832.  
833.  

834. //=============================================================================

835.  
836.  

837. /*

838. * The FLine class represents a line. Here, it is only used to call

839. * its intersect() method to compute the intersection of two lines.

840. */

841. class FLine

842. {

843. /* The line is under the form ax+by=c */

844. double a;

845. double b;

846. double c;

847.  

848. /*

849. * Creates a new FLine that goes through p and is supported

850. * by v.

851. */

852. public FLine (FPoint p, FVector v)

853. {

854. FVector vn = v.normal();

855. a = vn.x;

856. b = vn.y;

857.  

858. c = a * p.x + b * p.y;

859. }

860.  

861. /* Creates a new FLine that goes through both a and b. */

862. public FLine (FPoint a, FPoint b)

863. {

864. this (a, new FVector (a, b));

865. }

866.  

867. /*

868. * The methode computes the intersection point of two

869. * lines. We use the JAMA matrix package (see

870. * <a href="http://math.nist.gov/javanumerics/jama/">

871. * http://math.nist.gov/javanumerics/jama/</a>

872. * for details).

873. */

874. public FPoint intersect (FLine l)

875. {

876. double[][] tabA = {{ a, b},

877. {l.a, l.b}};

878.  

879. Matrix matA = new Matrix (tabA);

880.  

881. double[][] tabB = {{ c},

882. {l.c}};

883.  

884. Matrix matB = new Matrix (tabB);

885.  

886. double[][] tabX = (matA.solve (matB)).getArray();

887.  

888. return new FPoint (tabX[0][0], tabX[1][0]);

889. }

890. }

891.  
892.  

893. //=============================================================================

894.  
895.  

896. /*

897. * The MesMath class is an abstract class that is used to perform

898. * division related operations on integers. I haven't been able to tell

899. * java that I needed my divisions to be rounded toward minus infinity

900. * rather than zero, so I've had to do it myself.

901. */

902. abstract class MesMath

903. {

904. public static int div (int a, int b)

905. {

906. if (a >= 0)

907. return a / b;

908. else

909. return (a - b + 1) / b;

910. }

911.  

912. public static int mod (int a, int b)

913. {

914. return a - div (a, b) * b;

915. }

916. }

 

 

 

4.  Application examples

One of the main application of the Hausdorff distance is image matching, used for instance in image analysis, visual navigation of robots, computer-assisted surgery, etc. Basically, the Hausdorff metric will serve to check if a template image is present in a test image ;  the lower the distance value, the best the match. That method gives interesting results, even in presence of noise or occlusion (when the target is partially hidden). 

Say the small image below is our template, and the large one is the test image : 
 

We want to find if the small image is present, and where, in the large image. The first step is to extract the edges of both images, so to work with binary sets of points, lines or polygons : 
 

Edge extraction is usually done with one of the many edge detectors known in image processing, such as Canny edge detector, Laplacian, Sobel, etc. After applying Rucklidge's algorithm that minimizes Hausdorff distance between two images, the computer found a best match : 

 

 

For this example, at least 50 % of the template points had to lie within 1 pixel of a test image point, and vice versa. Some scaling and skew were also allowed, to prevent rejection due to a different viewing angle of the template in the test image (these images and results come from Michael Leventon's pages). Other algorithms might allow more complicated geometric transformations for registering the template on the test image. 

In spite of my interest for the topic, an online demo is definitely beyond the scope of this Web project !  So here are some Web resources about image matching with Hausdorff distance : 
 

 

• Hausdorff-based image comparison 
  A vision project that uses the generalized Hausdorff distance
  to locate a template image in an indoor scene.
   
• Hausdorff distance for object tracking and classification 
  Hausdorff distance used with a multi-layer perceptron 
  for recognizing and tracking cars in a road scene.
   
• Computer-vision-enabled ophthalmic augmented reality environment 
  An interesting application of Hausdorff distance for merging
  partially overlapping images into a single image.
   
• Comparing 2D images of the 3D world 
  Conference material by Daniel Huttenlocher on object recognition.
  Briefly introduces some problems of image analysis, and compares
  the generalized Hausdorff distance with more traditional methods.
   

---

Glossary

References