Hausdorff and Modified Hausdorff distance implemented using KDTree

Hausdorff KDTree.ipynb

{
 "metadata": {
  "name": "Hausdorff KDTree"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%matplotlib inline\nimport numpy as np\nfrom matplotlib import pyplot as plt",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from scipy.spatial.distance import cdist\n\n# these are the brute force methods\n\ndef hausdorff(A,B):\n    \"\"\"compute the Hausdorff distance between two point sets\"\"\"\n    D = cdist(A,B,'sqeuclidean')\n    fhd = np.max(np.min(D,axis=0))\n    rhd = np.max(np.min(D,axis=1))\n    return np.sqrt(max(fhd,rhd))\n    \ndef modified_hausdorff(A,B):\n    \"\"\"compute the 'modified' Hausdorff distance between two\n    point sets as described in\n    M. P. Dubuisson and A. K. Jain. A Modified Hausdorff distance\n    for object matching. In ICPR94, pages A:566-568, Jerusalem, Israel,\n    1994.\n    http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=576361\n    \"\"\"\n    D = cdist(A,B)\n    fhd = np.mean(np.min(D,axis=0))\n    rhd = np.mean(np.min(D,axis=1))\n    return max(fhd,rhd) ",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from sklearn.datasets import make_circles\n\npts, c = make_circles(noise=0.025)\nA = pts[c==0]\nB = pts[c==1]\n\nplt.gca().set_aspect('equal')\nplt.scatter(A[:,0],A[:,1],color='blue')\nplt.scatter(B[:,0],B[:,1],color='red')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 76,
       "text": "<matplotlib.collections.PathCollection at 0x7f5935527f90>"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEACAYAAACtefPrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4U2X2x7/p3nRhUSzYVgottIC0FgXcKh2hCigVf4hW\nXBCXcZRl3HDGZQRmsIILMyoOKCJWR8CNbWRRXIrIVmURtciihelCK5ZFutA07fn9cQxJmqSk994k\nN+n5PE+e9iZv3ntyk/fc9z3vWQxERBAEQWgjQb4WQBAE/0SUhyAIihDlIQiCIkR5CIKgCFEegiAo\nQpSHIAiKUK087rzzTsTFxaF///5OXy8sLESHDh2QmZmJzMxMzJw5U+0pBUHQASFqO5gwYQImT56M\n22+/3WWbIUOGYNWqVWpPJQiCjlA988jKykKnTp1abSN+aIIQeHjc5mEwGLB582ZkZGRg5MiRKC4u\n9vQpBUHwAqqXLWdiwIABKC0thdFoxNq1azF69Gjs27fP06cVBMHTkAaUlJTQ+eef71bbpKQkqq6u\ndng+OTmZAMhDHvLw8iM5OVnRuPf4sqWqquq0zaOoqAhEhM6dOzu0++mnn0BEunpMmzbN5zKITIEl\nlx5l+umnnxSNbdXLlptvvhkbNmzAr7/+isTERMyYMQONjY0AgHvvvRcffPAB5s2bh5CQEBiNRixd\nulTtKQVB0AGqlceSJUtafX3ixImYOHGi2tMIgqAzxMO0FbKzs30tggMik/voUS49yqQUAxGRr4UA\neEtXJ6IIQrtC6diTmYcgCIoQ5SEIgiJEeQiCoAhRHoIgKEKUhyAIihDlIQiCIkR5CIKgCFEegiAo\nQpSHIAiKEOUhCIIiRHkIgqAIUR6CIChClIcgCIoQ5SEIgiJEeQiCoAhRHoIgKEKUhyAIihDlIQiC\nIkR5CIKgCFEegiAoQpSHIAiKEOUhCIIiRHkIgqAIUR6CIChClIcgCIpQXatWEGwpLwc++wwwGoFr\nrgEiI30tkeApVM887rzzTsTFxaF///4u20yZMgW9evVCRkYGdu7cqfaUgk7ZsQPo0weYOBGYMAEY\nMACoqfG1VIKnUK08JkyYgHXr1rl8fc2aNThw4AD279+P1157Dffdd5/aUwo65Z57gJMnWWHU1AAl\nJcDcub6WSvAUqpVHVlYWOnXq5PL1VatWYfz48QCAwYMH4/jx46iqqlJ7WkEDPvsMmDIFmDYN+OUX\n9f1VVtofNzQAhw6p71fQJx43mJaXlyMxMfH0cUJCAsrKyjx9WuEMvP02kJsLvPwy8MwzQHo68Ouv\nju2IgGefBRISgPPOA+bNc93nkCFAeLj12GgEhg513nbVKuDii4GBA4HFi9V9FsE3eGW3hYjsjg0G\ngzdOK7TCo48CdXX8f2MjcPw4UFDg2O7f/wb+/nc2hJaWAo884nqwv/oqkJUFBAcDISHA1KnADTc4\ntvv4YyAvD9i2DfjmG17uLF2q3WcTvIPHd1vi4+NRWlp6+risrAzx8fFO206fPv30/9nZ2cjOzvaw\ndO2XU6fsjxsbnRs3334bqK21HtfVAf/5DzBunGPbmBhg/XperoSEsBJxxty5QH29fZ8vvcQKRfA8\nhYWFKCwsVN2Px5VHbm4u5s6di7y8PGzduhUdO3ZEXFyc07a2ykPwLDfeyIrBMogjIngZ05KYGPtj\ngwHo0KH1vm2XLs4IC3N8LjTU+v8PPwBlZUD//sC557bel9B2Wt6YZ8yYoawjUkleXh5169aNQkND\nKSEhgRYuXEjz58+n+fPnn24zceJESk5OpvT0dNq+fbvTfjQQpV1QVUX08cdEO3eq66ehgWjSJKL4\neKI+fbhPZ2zbRhQVRWQwEAUFEUVHE/3wg7pzb9lCZDQSsUWFKDLSev6HH+bjDh24zerV6s4lnBml\nY8/w+5t9jsFgcLCNCPZ8+SU7XgUH8zLjppuAhQt5NuBJvv+eZykhIcAddwC9eqnvc9s24J//BMxm\nYNIkIDsb+Ppr/muxxQBAVBRw4oTzJRAR8OmnvMszaBCQmqpervaI0rEnysOPiIuz31KNigI+/BC4\n+mrn7b/4gu0T0dHAn/8M9OzpHTmV8u67Vl8RC2FhQEUFcNZZ9m2bm4Hrr+ftZoMBaGpiQ+7o0d6V\nORBQOvbEPd1PaGoCjhyxf665Gfj5Z+ftly8Hbr2V7+JBQcCbbwI7d+pbgfTvz5/Tlg4dgM6dHduu\nWwd8/rm9Mff223mWIpt53kEC4/yE4GAe+LYDw2AALrjAefvHH7dO/5ubeSfllVc8L6ca+vblpUx4\nOM+qzjoLWLvWuTKoqODPZUttLS/nBO8gysOP+O9/eekSFcXT+WnTgEsucd62ocH+uLnZ/i6tV/74\nR55h7d7NCuLCC523GzSIbR4WgoLY5uFsJ0fwDKI8dE5BAdC9O9C1Ky89Skp4YP3yCzt6ueLOO9nD\n04LRyMsYfyAmhmdZrSmC9HR2SouIYENuSgqwZo33ZBTEYKpr1q0DxoyxLj+MRuDhh9nj80w0N7Pb\n+Ztv8vvy83mnJtBobubrEx3ta0n8F9ltCUAmTODBb0tKCrB/v0/EcY+mJmD7dvY+u+giXmN5ifp6\n4E9/Aj76CIiNZU/WQFSYWiO7LQFIx45sKLXdgWjp8ekTGhvtXUItNDQAw4YBu3axESIqCti8GUhK\n8opYd98NLFvGrvdHjwJjxwKbNgGZmV45fbtDbB465qGH+A4aEsI7DpGRwPPPe+nk770HjBzJnmjf\nfcfPbdnCxpfwcCAxkfd+bXnpJZ511NQAv/3Ghpm77/aSwBypaxuzYzLxbo3gGWTm4SOIeKy98AIr\nhqlT2dPSlsREHrdvvMFT8htu4OxcHmfBAuCBB9iYYDAAq1ezU8Xw4awUAA4+GTaM/1pyDf7wg33E\nW1MTsG8f/0/EYbn19UByMmtEjYmKsg/uCw1l5dtWTp3i6x4RAfTrx5MowQkq3eI1Q0eieIXXX7eP\n7zAaiQoKvCzE8eNEmzYR7d9v/3yPHlbBAA5suflmDjixfT4mhmj3buv75s2z/1ChoUTXXUdkNhON\nGUMUEcGBMqmpHKSjMYsXW08fHk503nn8EdtCWRlR9+780aKiiIYNIzKZNBdVVygde7oZse1NeWRl\n2Y9DgOjKK70oQFERK4PYWI5EmzzZ+lr37o7CTZjA7WyfCw8nqqy0vq+piSgvj583Gon69SM6coTo\nlVfslUpICFFurkc+1saNRI8/TvT880THjrX9/cOHEwUH2wftzZmjvZx6QunYk2WLj3A2nVYyxVbM\n//0f+3JbeOMNjskfNoxzE/7tb/Z7xH/+M1twX3vN6p316KPstWYhKAhYsgQ4fJiXJ927s8X366/t\no93MZkd7iUZcfjk/lFJcbG+grq9n+6/giCgPHzFjBlBYaD8+p03z0smbmzk1mC1NTcDevaw8HnyQ\nBSooYAeKf/wDyMgA5sxhBbNvHxsDLrvMef/dutkf9+/PdhGLPSQ4WLchsOnp7NlqNvOx0ejay7W9\nI34ePqS4GFi0iG/YEyYAaWlePHmPHsDBg9bjqCh2kPBE9jaTiUN/v/6aFUdMDG/hNjcDt90G/Pgj\nB7a89RbPVnxIZSVwxRU8eWpqYl26bJlH7Lu6QZzEhLbx7becndhk4sdDD7EbqqdobuZz1tdzNF9Q\nEHu8VVbyKA0O5rRh+/c7T0VWVMRhs4cPc2DLO+8A55zjEVFPnAA2bODdrgsuCPwoXVEegj2WAdka\ndXU8WLt08X6+v+3bgT/8wT55R0wMZzxqGSp8+DAvcyxtQ0J4KbRjh+ZibdwIXHstm3UaG9lL9a67\nND+NrlA69mQH28u0DCNXRUUF+1bYhtBWVQGDB3NUWUwMGzBdYTSyLcMXiUKjoqyGBQtms/MglU2b\nHNt9/z37n/Towe8ZNsx57Yg20NgIjBrFriwnT7K/x+TJwIEDqroNWER5eInFizmxTWgo2xlbJvZp\nE0S8I9KzJ8fkJyVZA15Gj+Y7siWJx113eWxnQxWpqex0Zgn9NRr5lp+c7Ng2NtY+/h7gz3fjjWy3\nqa3lGcu116oSqbLSMR9IaCjbpgQnaLBNrAk6EkVzvvnG3kUiNJT9PBTz3/+yB5OtE1d6OlFzM2cp\ntvXFiIggeuklzT6LppjN7C03ZQrRwoXsJ+KMxkaiSy+1+opERRFdfz1nY7b9rMHBRPX1isU5dcre\nHcXi5/H994q79AuUjr0AtiHrh40b7ZcrjY0cJuIWq1ZxnEnnzuxXkZDAvtO2QRxEvM1qqYtw7Jj1\ntZAQjxkWVRMc7J5BISSEE7K++Sa7uF98MS9d1q937E9FNqDwcJ4hjhvHMw6TiTOy9eunuMvARmMl\nphgdiaI577xjP1EAiM46y403zp9vvRUGBxN17kxUUUH0wQeOHaam8ns+/JBvl5GR3Obyy/nOHWiY\nTESDBvH1CQrivy++aH29rIzossv4+ZQUoq1b3e66ooLo888dvfYDFaVjT3ZbPIDZzLue69ez28LM\nmVwN7YcfrN6LS5YA1113ho66dmUDqIXQUO5s6lTetvzwQ77TBgfznTk9ndt9/z3w1Ve8i3LddYHr\npGAycU2Iigp2K/3DH/h5IvYb2b/fesFjYti5rWtX38mrUxSPPQ0VmCp0JIpqbrvNfsIQF0f0yy9E\nS5dymMd337nZUadOjmv6v/+dX2tu5upLGzcSnTjhsc/il1RWcnyN7bWLjSVasUJ1101NRA88wJO6\n6GiiJ57gr8KfUTr2ZOahMSYTbxzYxkdERwOvv86pMdrEAw9weLytD/u2bcD552smb0BSV8dxOLZb\nJ9HR7EHb0MCVorp142zLbcx0Nns2p4G0/Uqeew64/34N5fcy4uehc9r03RQUAPHxPCXv25c9MQcO\n5PLyojjOjNEIPPEEKwZLRrOLL+YIt+uv59H++OOcJtE2/4gbLF9uH+NXV8fu6+2RAF0M+46wMLbW\nf/gh/7CCg/m3O3y4mx188gnfxiy/0FOn+Pi55zwmc0AybRo7y339NRuebr6Zd6IsyuLUKd65WbYM\nuOUWt7vt0oU3tSw3g6Ag+8Di9oQsWzyAxWD6ySf8u332WZ5IuMUf/8hLFVuSkrjmgqCc5mbW7Lbr\nychIrjJ1771ud/Pjj6yTGhpYiUREsE9ejx4ekNlLSAJkHRESAjz1FD/aTKdOOs167OcEBXEg4IYN\nVnd+y3NtIC2NN7OWL2flccMNjhkI2guqbR7r1q1DWloaevXqhdmzZzu8XlhYiA4dOiAzMxOZmZmY\nOXOm2lMGNg88wAokNJR/nUYj3x0F9bz/PgevdO4M9O7NVaJSUtrcTWIiRwdMntx+FQegctnS1NSE\n1NRUfPrpp4iPj8fAgQOxZMkS9OnT53SbwsJCzJkzB6tWrWpdED9atjQ0cAxWXJwKF4oDB4CFC3lH\n4NZb7SNJKyvZaFpbywY+qR2ga6qr2VF22zZepr7xBtu5/QWfLFuKioqQkpKCpN/rcuTl5WHlypV2\nygOA3ygFd1ixgu1rROzOvHo1cOmlbexk717ePamt5bX4vHm8k2LJn9e1K/CXv2guu6A9RJznaPdu\nvg9UVfHXuH8/F+oOZFQtW8rLy5GYmHj6OCEhAeUt0tsZDAZs3rwZGRkZGDlyJIr9OESxvJwVR10d\nG+2PH+fSJrZhJm6Rn88Rr5aAl7o63joUfMNvv/GXqYBff2UbiMWlhIjNVZs3ayifTlE18zC4kWJp\nwIABKC0thdFoxNq1azF69Gjss9TyaMH06dNP/5+dnY1sT6TEU0FxsWOhtKYm3vHr1asNHf32m6Pj\nh21SHME7mM18N1i+nI8tOQcjItzuIiLCMUdLc7NXq2y2mcLCQhQWFqrvSI1b65YtW+jqq68+fZyf\nn0+zZs1q9T1JSUlUXV3t8LxKUbzCjz86rz7QZu/wZcvsY7+jogI/v78e+cc/7L+HiAiiBx9sczdT\npli7iYggGjjQv2IRlY49VSO2sbGRevbsSSUlJdTQ0EAZGRlUXFxs16ayspKaf3f+37ZtG3Xv3t25\nIH6gPIiInnqKfyiWcicLFyrsaOFCoqQkooQEovx8/w+Q8EeGDXOsTzNgQJu7aW7myOn77yd64QVV\nKUV8gk+UBxHRmjVrqHfv3pScnEz5+flERDR//nyaP38+ERHNnTuX+vXrRxkZGXTJJZfQli1bnAvi\nJ8qDiIukrVjRfkK2A5b77uPMTLbFqMaO9bVUXkfp2BMPU6H9Ul3Nu16W3KfR0ezO7rY7cGAg2dMF\nQQl1dVx9q7kZGDKkXXrzivLwF8xmvuOdffaZSyMIgheQkHx/4OOP2fU8KYnDM1uWFBD8GrOZk5qZ\nTL6WxDuI8vAWR44AY8awc9ipU5ykeORI9jIVfMvGjezjcemlnGRZAdu2cZ7plBTOQ2RxHQlkJKrW\nW+zZ4xgIQ8Sh9pLgx3cUFbF/uSXPx7ffsrvoPfe43UVDA+drsXVSvfVWjkJISNBYXh0hMw9vkZBg\nX9kN4Pltew7L1AOvv26fTayuDpgzp01dlJU5Lxb1ww8ayKdjRHko4Phxjub+/HPHH41LevYEHnuM\nQ+xjYzkRzfPPB370lN5xFmIR1LZhERfnWDnTZALOO0+FXH6A7La0kQMHuMKjycS7e716cZUDS9VE\nB9av51Dcs84CJk3ios3793PMtj/FbQcq337L9T8ttiejEXjtNWtqwuXLOeOx2czf3x//6FThvPkm\nZ4sMDeWmDz7IVTL8ASm94EH27iX6+GOi0lKi7Gz7io4REURPP+3ijYsWWYMeQkK4BsORI94UXXCH\nb74hGjOGaORIouXLrc+vW2cfzGQ0Ei1Y4LKbffs4bGnXLi/IrCFKx57MPM7A00/zIyyMZxtRUY7F\n2MePd2Gkb1m0KTwceOYZvi0J+mfMGMfU6BdeCHzzjW/k8RCSw9QD/PgjK476eqtNzWSyKhKAZ7mX\nXeaig5aJPsxm2Zr1J5yF5oeHe18OnSIG01YoKXGsmxweDvTpw39DQzmjv8tazTfdxIZR2zePGuUx\neQWNefRR+8QcRiNgk3PGXXQ4odYEWba0wsGDbNO03cnr0IFTjJ44wTemDh1a6cBkAh55hKe+HTsC\nL77Y5mzdgo/ZvRt46SWr74clVaQbHDwI5Obylm3nzlyfeNgwz4mqFIlt8RDvvAPcfTf7dwUFccXC\nrCw33/zLL5zkNDiYZxydOnlUVkEDqqv5znDeeaoKhBPxTlxJiTXTWFQUK5Lu3TWSVSNEeXiQmhre\nYU1MbEOGup9+AgYNsjqGRUUBO3cC557rMTkFlfz1r1zmIjSUt9YLCxVXc6qu5q/aNs4lNpYT5t9w\ngzbiaoUExnmQ6Gi+i7QhtSXw8MPsTVZby4+jR4Enn/SYjIJK1q4F5s7l0V5by26jY8cq7s5ZZH9z\nM8dDBgqiPDxFebl9ZlyzmX+Qgj7Ztct+d6y5WZV/eVgYe7kbjXzTiY5mc9cVV2ggq06QrVpPMXw4\np1u3FKw2GoERI3wrk+Canj15lNtupduUFVHCxIlWt5CEBDaeulFwwG8Qm4cLjh7lSvcmE9s62xyn\nYDazdf4//+HjP/2Jd1vaGDcheInmZl6mfPyx1VD6+efAgAG+lcsLiMFUQw4f5uqPlrpMoaEcv5Ke\nrqCzpia+3YjS0D9EXPL+6FGeMnTurEm3tbUcSNnQAOTkcCCdnhDloSGTJgGvvmqNlDQYgOxsvhEJ\nQls4doz10JEjfBwSAmzZAqSl+VYuW2S3RUMqKuxDrInYZUMQ2sozz7DtvKaGHydOcPRtICDKwwm5\nufZeyZGRwLXX+k4eQWd89RVrgEceYS+wVjh0yN7XgyhwNt1EeThh/Hj+XRiNvOWWlwf84x9neFN1\nNeeGOHHCKzIKPuKjjzht4bx57FB2wQXAzz+7bJ6TY38jiojQp4u6EsTmoQVvvQXcey9rmqYm4P33\nZVs2UElPB777znocFARMngz8619OmxOxv+DLL/P/I0YA771nHy/pa8Rg6itKS4HUVPvouagojp6L\njvadXIJnSEnh0ANb7rqLc6G2gtnM9xU9RvSLwdRX7N/vGLdvMLBSEQKPO+6wzzkZGcmp0s9ASIg+\nFYcaxMNULT17Olb5MZk4TmLgQOD228XHI5B4/HFefyxaxAaM/Hzex2+PKEpeaMPatWspNTWVUlJS\naNasWU7bTJ48mVJSUig9PZ127NjhtI0GoviOV17hZKYdOhAFBxOFhXHOy6goohtvJGpu9rWEgg9p\navK1BK2jdOypGrFms5mSk5OppKSETCYTZWRkUHFxsV2b1atX04gRI4iIaOvWrTR48GDngviz8iAi\nKisjWryYlYglYS7ACXT37/e1dIJW1NYS3XEH0bnnEmVkEG3e7LLpt98SJSURGQxE8fFE27Z5Uc42\noHTsqZpPFxUVISUlBUlJSQgNDUVeXh5Wrlxp12bVqlUYP348AGDw4ME4fvw4qmyTAusQImD2bM4S\nFhXFYSkt63I4EB8P9O7Nvuy2hISwd5AQGNxyC7B0KXsSfvst78W2NKCC7edXXsnZxIjYUeyqq+yr\nyvk7qpRHeXk5Em0iDxMSElBeXn7GNmU695JZsoRLdfz2GwfFvv028NRTbryxb19O5GCxcQQHswbS\nky+yoBwi9vNoGbr/yScOTffvd17wOpCqyKkymBrcjC+mFttArt433Sa5bHZ2NrJ9ZIhascIaSQ/w\n/ytWsG2sVSIj2fvwtts49XqfPqx52pRFSNA1lqpOFoKCnFb86tLFuR39nHM8LJ8bFBYWorCwUHU/\nqpRHfHw8Sm22JEtLS5HQorJvyzZlZWWIj4932t90BZmpPUFcHK82bH8jZ5/t5pt79GAFIgQeBgO7\nGj/1FN9RwsJYG4wZ49C0Wzf2Uv7Xv9i/IzgYmDCBM9L5mpY35hkzZijrSI2hpbGxkXr27EklJSXU\n0NBwRoPpli1b/MJgWl5O1KUL2z7Dwoiio4m2b/e1VIJuWLmS6P77iWbOJDp+vNWmX3xBNHcu0fr1\n3hFNCUrHnmoP07Vr1+KBBx5AU1MT7rrrLjz22GN49dVXAQD33nsvAGDSpElYt24doqKisGjRIgxw\nkmBFbx6mR46wG7HJBFx3HbtzuM1HHwFffw0kJfESRkUWbsGHEPEPINC8u1og7ul64cknea5aV8dr\n4Usu4exU4ijmX7z+OsesmEycTeyjj/SXxUcjRHnogZoazj7V2Gh9Ljqaa7cEUubbQGfrVs5WbLGa\nh4QAF18MbNzodhdNTcCsWZzK8uyzgeeeAzIyPCSvSqRWrR6oqWHLmK3yCAqSMH1/46uv7L9DsxnY\ntq1NXTzyCPDaa1b9c/nlXHxOYRkYXSJz6Vb47bc2+nfFxbGdIzjY/vnBg7UUS/A0Xbs6Bju2MZ/p\nwoX22/0mE1cdDSREeTjh1Cngmmu4aFinTsC4cW54mAK8lffZZ8Bll7GzWFoaJz7Vw+a+4D433cRJ\nfqKj2cXYaATefLNNXbQ0cRkMjvcUf0dsHk546CFg/nxrig6jEfjb37gaodBOMJvZVlVdzWuO3r3b\n9PaZMzl/aV0dK5IOHdi7tFs3D8mrAjGYasiAAVxW1parruJNk1ZpbgZmzOA5a1gY/4LGjfOYnIJ+\nIQIKCoAPPuCJ57Rp+itwbUEMphqSnMzGraYmPg4L4+fOSH4+8Pzz1sXuPffw2ufqqz0mq6BPDAbO\nG3THHb6WxHOIzcMJ//wnxybExPAjPt6NBMgAx7G0DIpZvNhjcgpeYPdu4KKL2Ig6ZkxghcWqRGYe\nTkhIAPbuBQoLeb165ZVOY58caZmzNCgIiI31hIiCN6iqYv8cy1b7Rx+xJX3TJt/KpRPE5qEln37K\nvux1dWxaj45m40kgbe63J95/n5MbnzxpfS4khI2oAXRTkATIemDYMJ6uPPQQb818+60oDn/GaGTL\nZ0vOEOvy3Xf8U8jI4F06t7b5/RBZtrTCgQOcAapvX7aBuMXAgfwQ/J+cHC618OOP7PwTFcWFjFtR\nHocOsZtPTQ3rnQMHOMhy/nwvyu0lZNnigieeYMNpWBjfOVasCJxKX0IbqK8H/v1vLit5xRXA2LG8\nleKCl18Gpk4FGhqsz0VG2tvR9YZs1WrI119zYGx9vdVRbMwYrniuKjiWiLdtXnuN187TpnGGGEFf\nlJTwD+DkSfbTefhht98aEuL4Gwk0z9LTKMoC4gF0JAotXkwUE2OfBD0sjOjoUZUdz55NZDRaOzUa\niVat0kRmQSMOHuQSGkFB1u9oyRK33/7LL0Rnn80VOCxvnz7dg/JqgNKxJwZTJ/Tr52jkiokBOnZU\n2bEzP5C331bZqaApr77KBovmZj6uq+McLW7SpQuwaxdv0uTmcu0vt5Jn+yGybHFCejrHJfzlL2zz\nCA7mMAc38z27pqUfiMHAQQ+Cfqivt7oWW7A1YLhBfDzroEBHDKatUF3NfkI9evBMpKICOO88FRXO\nN2wARo60RktFRwPbt7NFX9AHRUXAH/5gnSEajWzz+PvffSuXB5HAOA+yeDFw9908AzEYgJUr+fel\niB07uDBMeDjPbcUPRH988glPO2tqOAftk08GdBpJUR4eorQUSE217roAbP+oqlIwAzGZuBTdtm3A\n+eezB1FUlKbyCr6lpASYM4c3am691T+290V5eIj163lr3zaTYHQ0TyDaVIODiJcsGzawJoqI4KJQ\nRUWSXT1AOHiQcwidPMn21shIYNEizi2kZ8Q93UP06OFY+aupSUFSl7Iydl23TGFOneKahNu3ayGm\noIaDBznBcUwM0L8/8P33irp59VWr4gD4q27DRo3fIcrjDKSkANOn812kQwfr3aTlxskZMZud56YL\n1MAHf6GxERgyhD0Da2o43deQIYqSVtfVWRWHhTZu1PgVojzc4NFHOdjpgw+AffsUTkO7d+dliiUu\nIjSUc/JfeKGmsgpt5OBB3lazjHoinlru2tXmrm6+2T51g9EI3HmnNmLqEbF5KMRsZoVCxDPd0FA3\n3nTiBPDgg3yX69sXeOmlgC0k5DdUVbFit50iGI2cs+OCC9rc3ccfc0B1bS1v1DzxhP43asRg6kVq\nanhmu298RaVZAAAUIElEQVQfH3fvzqU+VHugCr5h8mRei1qq/F11FVdrUu0V6B+I8vAiDzzAIdaW\nm1VYGN9lXn9doxNs3QoUF/Me8WWXadSp4BIiDpvetYu30MaN0/90QUNEeXiRK68EvvjC/rmBA3nX\nVTUzZgDPPst3PSK+K86apUHHAn74gbODhYYC48dzvknB+8rj6NGjuOmmm3Do0CEkJSXhvffeQ0cn\n8/akpCTExsYiODgYoaGhKHIxwvxJeUydygFPp07xcXg4R9bPm6ey44oKTtNu6Rhgf5DiYvFEVcuW\nLeyxdeqUNTRgxw65rvCBn8esWbOQk5ODffv2YejQoZjl4u5oMBhQWFiInTt3ulQc/saMGVzbxWhk\nB9F+/dhxVDW//OJY5jA8HKis1KDzds7Uqda9VLOZa4lqOKMjAt54g+tjjx3Lk5yAR1EgPxGlpqZS\nZWUlEREdPnyYUlNTnbZLSkqiX3/99Yz9qRDFJzQ1Ee3ZQ1RcTGQ2a9TpyZNEHTvaJxKJjSU6dkyj\nE7RjUlLsrytANHasZt0/+6w1VYvBQBQdTXTggGbdexSlY0/xzKOqqgpxv28zxsXFoaqqymk7g8GA\nYcOG4aKLLsKCBQuUnk53BAVxKdo+fTTMFBUdzUFZXbvyCbp0AdaulW0ctaxfz8lFbYmM1NRvfM4c\nayAuEXuXvvOOZt3rklaDKnJyclDpZMr89NNP2x0bDAYYXGxrbdq0Cd26dcORI0eQk5ODtLQ0ZGVl\nOW07ffr00/9nZ2cjOzv7DOIHIGedxWkKu3XjYkOCeiZOZE9SW3JyOLekRjgzGejVhFdYWIjCwkL1\nHSmd6qSmptLhw4eJiKiiosLlssWW6dOn0/PPP+/0NRWiBA5vv00UGclLlchIoief9LVEgUGXLo5L\nlocf1vQUTz9tn2EyKopo3z5NT+ExlI49xcuW3NxcFBQUAAAKCgowevRohzZ1dXU4+XvBnNraWnzy\nySfo37+/0lMGNrW1XNu2vp6NefX1wAsvtBPLm4e59lretbJgNHKEs4Y89hjw3HPApZdai8q1Kera\nH1Gqraqrq2no0KHUq1cvysnJoWO/G/XKy8tp5MiRRET0008/UUZGBmVkZFC/fv0oPz/fZX8qRAkM\nfv6Zb1e2d8cOHYhWr/a1ZP5JYyMbmpubierqiMaN46nBWWcRvfGGr6XTFUrHnjiJ6QWTieNcbAsp\nG43s49G9u+/k8kcKCoB77+Vt2YQENpgmJ6vulohNJy130/0dyefhY3bt4ppAaWmc8rJlDpAzEhbG\nOyudOrHSiIwE3nyTq5UNHMjRdy++qF8rnF7YvRu47z6OHWhs5KjZESNUd/vf/1pTMqSlAT//rF5U\nv0fD2Y8qdCRKmzl0iPf1LauNyEii225T2JnJRPS//xHV1xN99RV3ZmuFc2FwFn7n9dcdl39BQXw9\nFXLggL0x1GAg6tmTV0SBgNKxJzMPDfjoI/ucPvX1wNKlCicJoaFAYiIb+AoK7JOn1ta2j5z+anAW\nr2I0nrE4dWt8/bW9Lw8R57ZVkC8ooBDloQFhYY5BmCEhGkR0h4c7dhJoC26tueoqXqZERXFaQaOR\nvbVUfBlduzpmCAsK4u7bM2Iw1YBjxzi3z6+/8gzEaOSEMH/7m8qO9+5lR7HaWr7dGY3AW285OjeV\nlvIAMZuBG28EevdWeWI/4csvOZQ5MRG44Qbr9ICI88VWVrK9SGVdHCKOV1m3zroW+ve/gTvuUP0J\ndIGE5PuYqiqOs6qoAEaNAm65RaNcMnv2sO9zbS2H7ubk2L9+4IC9gomI4EE1YIAGJ9cxL7zAdRwt\n2x9DhrBV00N5OIg4S1h5OV/ujAyPnMYniPJor9x+O886bOfVQ4cCn37qO5k8zalTQGysvct5dDQr\nD41CGsrKgJdf5mzoN96oWbe6ROnYk4Ih/s7Ro44L8mPHnLetqWE7ilsJV3VMTY3jDCMoiBMZa0BZ\nGdcr/u03zoVcUMC75mPHatJ9wCAGU3/nppscU3bn5dm3OX6c0xlafEj8vWz7WWex45ztFkhzMzB4\nsCbdv/aaVXEAHC372GOadB1QiPLwd267jYswd+kCdO4M/PnP7KVmy4QJwDffsEHVbGYbyrJlvpHX\nFRs2cOGlfv2AZ55xnE3ZYjDwsiwzk2dRCQnA6tWapRWsqbEqDgu2O+YCIzYPL1JWxjWum5t5c0AD\nj2n3OOcc4MgR++eGDuW6AIordmvIzp3A5ZfbV6afOpWrbfmALVv48lgUhtEITJqkUbY4HaJ47Kl0\nTtMMHYniEQ4c4Di30FCikBD2SN21y0snz8hwDEkPC2NPzPvv95IQrfDII47ynXuuY7vFi4kyM4kG\nDCB6911NTn3yJNH773PXR45Yn1+9mqhvX6Lu3Ykee0zDbHE6ROnY082IDXTlcdtt7CVtOz66dvVS\nzoft24liYux96C2P4GCONE1PJ9q2zQvCOOGJJ1gOW7m6d7dv89579j7iRiPRihWqTnvkCNF55/Fl\niY7my+AvqQO1RJSHzhk+3HHcWlKU/u9/XhCgooJozhyiiAjnggA8gg4d8oIwLSgp4Qth0a5GI9Gi\nRfZthgxxlHfYMFWnnTSJZ4K2ITC/Z5NoVygde2Iw9RI33cQe0y05dQp47z0vCNCtG3D//a1X6CYC\nPvvMC8K0ICmJyyDccw9fqHffdXTfdOaWb5vgRwGHDtm7ijQ3s7Ou4B6iPLzE+PHAk086ep0Stb6x\noCnh4cDnn7M7tzMMBtfKhQh45RW2JI4bp31MenIyl+FbupQzf7XkyScdt6T/+ldVpxw2zL7LyEj+\neIKbaDwDUoyORPEojz1mn6I/JoZn7a3R1ET0228aC9LQYJ94MyKCKC3Ndei6reBBQVwioqLC+vrH\nHxNNm8Yh8SaT8z7UWh03beKMYOPGEW3dqq4v4ut6331sbgkOJho9mujUKdXd+h1Kx55uRmx7UR7N\nzUQvvEB0ySVE115L9MMPrbdfupRTeoSEECUlEe3dq7FAK1cSTZnChUdqaly3a2lsDQ8neuklfu2Z\nZ1ixGAz897LL7BVFeTnvkAQFsW3j/fc1/hCuMZnOnHfDZGqfSsOCKI8AZM8exyQ0SUk+SkLTMsFO\neDjRv/7FI8/W6mgxvK5da31vZqb9bkpkJNF333lU3KoqokGDWF9FRvKESHCO0rEnNg8ds327YxKa\nsjIO1vI6f/qT1UBgMLD9ZMwY9qRq6WBkMFgz5ZjNnKPR1mXTYAA2b/aouGPHsg22uZlFnDKFnb8E\n7RDloWPOPdfRmBoa2vqGicd49ln2+Lz4Ys45sG0bu4PHxgLnn8/ZjywQcSwNwNqv5TZTUBAne/Yg\n27bZZ3drbAS++sqjp2x3iPLQMdnZQG6uNSmWJSeyh1JWtE5QELuMb9kCrFzJWYAtrFsHZGWxVktO\n5sQXljgTgwFYsICFj4zkNoMGOd9R0ZBOneyPw8I4I5igHRLbonOIOGasvBy48EL7MetX7N7NS5W4\nONaICgv8NjcDv/xizWTuinXreFVlMPCjf3++jv6ejcATSDKgdg6RRpnLdExJCfthHD7MJpRZs4CH\nHnLdft8+YONGDjYeNcp+ZSVYEeXRTlm1iiPuT5zg9HgrVgTu9Dw9natvWuxARiOvkC6/3Ldy+TtS\n9KkdsmcPcPPNnEysqYl3Z0aN8rVUnoHIXnEA1s8s+AZRHn7Mxo32x2Yzb0/axmtY2LOH68scOOAd\n2bTGYOC0JLaEhADnnecbeQRRHn5Nly6OOy8REY5r+2efZWPrrbfy1H/BgjP3vXs3z2KysjjkRA8r\nyiVLeOcpNpb/5uQA113na6naMUq90t577z3q27cvBQUF0fbt2122W7t2LaWmplJKSgrNmjXLZTsV\norRbGhuJsrLY+TM8nD0p33rLvk1JiWMUfkQE0a+/uu533z57b/SoKKL8fK5++cEHZ47F8SRlZUTL\nl7MsgVLu0dcoHXuKR+yePXto7969lJ2d7VJ5mM1mSk5OppKSEjKZTJSRkUHFxcXOBRHloYjGRo5/\nefFFom++cXx9wwbOYGarPGJiiHbvdt3ntGmOiYvCw1mJxMayy/zq1cpl3r+flVBRkfI+BO1QOvYU\nb16lueFwUFRUhJSUFCQlJQEA8vLysHLlSvTp00fpaYUWhIRwCgxXpKY6t4H06OH6PRbfCFtMJi48\nb2HcOK7w0Nbt4aVLgbvuYrnNZt4pmju3bX0I+sCjNo/y8nIk2uSOSEhIQHl5uSdPKbQgLo4HbFQU\nb23GxvL2bmsu7rfcwg5YFsUQFuboXHXyJCsUW9as4RQbc+dykqOWmEysLOrquLRBXR2waBEXknbG\nhg3AzJlso7FVXII+aHXmkZOTg8rKSofn8/PzMcqNPUFDoHst+QmjRnE9pKoq9gE5U63slBT2Qn/q\nKS75kpUFPPec9XWDAejZ077wfH4+8PTTrBAiI4E33gC2brU/l7NaVCEhwP/+xyVlbZk3D3jkEVZC\nERFstN2yRep864lWlcf69etVdR4fH49Sm7xupaWlSGiltsZ0m1T72dnZyA7kGn9eJjy8bdua559v\nX9olPh6YPJkVR1wczzIsNDVxzJxleVRfD+zfzy7iubnWdl26sFu57azEbAYuuMD+3ERcesZS+qCu\njr1FV63ikhWCOgoLC1FYWKi+I7XGluzsbPrGmaWOiBobG6lnz55UUlJCDQ0NYjD1cxoaOE9Gy12O\nujrH5OfR0URvv+3Yx7ffEp1zDu/4REQQLVni2Kax0dFgazQSLVjgmc/V3lE69hSP2GXLllFCQgJF\nRERQXFwcDR8+nIiIysvLaaRNCuo1a9ZQ7969KTk5mfLz810LIsrDr7nsMi4FY6s8XGWFb2oiqqx0\nna2QiLegbXMMGY1eKlPRDlE69iS2RdCEo0c5yfOmTWxXWbRIXenYo0d5R+err7g07RtvSHJiTyGB\ncYIgKEIC4wRB8CqiPARBUIQoD0EQFCHKQxAERYjyEARBEaI8BEFQhCgPQRAUIcpDEARFiPIQBEER\nojwEQVCEKA9BEBQhykMQBEWI8hAEQRGiPARBUIQoD0EQFCHKQxAERYjyEARBEaI8BEFQhCgPQRAU\nIcpDEARFiPIQBEERojwEQVCEKA9BEBQhykMQBEWI8hAEQRGiPARBUIQoD0EQFKFYebz//vvo168f\ngoODsWPHDpftkpKSkJ6ejszMTAwaNEjp6QRB0BmKlUf//v2xfPlyXHHFFa22MxgMKCwsxM6dO1FU\nVKT0dD6hsLDQ1yI4IDK5jx7l0qNMSlGsPNLS0tC7d2+32iqpwK0H9PhFi0zuo0e59CiTUjxu8zAY\nDBg2bBguuugiLFiwwNOnEwTBS4S09mJOTg4qKysdns/Pz8eoUaPcOsGmTZvQrVs3HDlyBDk5OUhL\nS0NWVpYyaQVB0A+kkuzsbNq+fbtbbadPn07PP/+809eSk5MJgDzkIQ8vP5KTkxWN/VZnHu5CLmwa\ndXV1aGpqQkxMDGpra/HJJ59g2rRpTtseOHBAC1EEQfASim0ey5cvR2JiIrZu3YprrrkGI0aMAABU\nVFTgmmuuAQBUVlYiKysLF1xwAQYPHoxrr70WV111lTaSC4LgUwzkatogCILQCj7zMNWjk5m7Mq1b\ntw5paWno1asXZs+e7VGZjh49ipycHPTu3RtXXXUVjh8/7rSdN66TO597ypQp6NWrFzIyMrBz506P\nyNEWmQoLC9GhQwdkZmYiMzMTM2fO9LhMd955J+Li4tC/f3+Xbbx9nc4kk6LrpMhSogF79uyhvXv3\nntHgmpSURNXV1bqRyWw2U3JyMpWUlJDJZKKMjAwqLi72mExTp06l2bNnExHRrFmz6C9/+YvTdp6+\nTu587tWrV9OIESOIiGjr1q00ePBgj8njrkxffPEFjRo1yqNytOTLL7+kHTt20Pnnn+/0dW9fJ3dk\nUnKdfDbz0KOTmTsyFRUVISUlBUlJSQgNDUVeXh5WrlzpMZlWrVqF8ePHAwDGjx+PFStWuGzryevk\nzue2lXXw4ME4fvw4qqqqfCoT4L3fj4WsrCx06tTJ5evevk7uyAS0/TrpPjBOb05m5eXlSExMPH2c\nkJCA8vJyj52vqqoKcXFxAIC4uDiXPzJPXyd3PrezNmVlZZrL0haZDAYDNm/ejIyMDIwcORLFxcUe\nk8ddvH2d3EHJddJkq9YVenQyUyuTwWBQfO62yvT00087nNvV+T3tjOfu52559/LE9WpL3wMGDEBp\naSmMRiPWrl2L0aNHY9++fR6TyV28eZ3cQcl18qjyWL9+veo+unXrBgDo0qULrr/+ehQVFakaFGpl\nio+PR2lp6enj0tJSJCQkqOqzNZni4uJQWVmJrl274vDhwzjnnHOcttP6OrXEnc/dsk1ZWRni4+M1\nk0GJTDExMaf/HzFiBO6//34cPXoUnTt39phcZ8Lb18kdlFwnXSxbXK216urqcPLkSQA47WTWmgXb\nGzJddNFF2L9/Pw4ePAiTyYR3330Xubm5HpMjNzcXBQUFAICCggKMHj3aoY03rpM7nzs3NxdvvfUW\nAGDr1q3o2LHj6SWXJ3BHpqqqqtPfZVFREYjIp4oD8P51cgdF10mp9VYty5Yto4SEBIqIiKC4uDga\nPnw4ERGVl5fTyJEjiYjop59+ooyMDMrIyKB+/fpRfn6+z2UiIlqzZg317t2bkpOTPS5TdXU1DR06\nlHr16kU5OTl07NgxB5m8dZ2cfe758+fT/PnzT7eZOHEiJScnU3p6utthC56Uae7cudSvXz/KyMig\nSy65hLZs2eJxmfLy8qhbt24UGhpKCQkJtHDhQp9fpzPJpOQ6iZOYIAiK0MWyRRAE/0OUhyAIihDl\nIQiCIkR5CIKgCFEegiAoQpSHIAiKEOUhCIIiRHkIgqCI/weqzWIO4C4eZwAAAABJRU5ErkJggg==\n",
       "text": "<matplotlib.figure.Figure at 0x7f59357e4dd0>"
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "hausdorff(A,B)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 77,
       "text": "0.28356701351218622"
      }
     ],
     "prompt_number": 77
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from scipy.spatial import KDTree\n\n# now here are KDTree based implementations\n\n# Hausdorff distance\nd1, _ = KDTree(A).query(B,k=1)\nd2, _ = KDTree(B).query(A,k=1)\nmax(np.max(d1), np.max(d2))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 78,
       "text": "0.28356701351218622"
      }
     ],
     "prompt_number": 78
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "modified_hausdorff(A,B)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 79,
       "text": "0.20149808169735792"
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Modified Hausdorff distance\nmax(np.mean(d1), np.mean(d2))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 80,
       "text": "0.20149808169735792"
      }
     ],
     "prompt_number": 80
    }
   ],
   "metadata": {}
  }
 ]
}

hausdorff_kdtree.m

n = 50;

theta = linspace(0,2*pi-2/n*pi,n);
noise = 0.025;
inner_scale = 0.75;

unit_circle = vertcat(cos(theta), sin(theta))';

A = unit_circle + randn(n,2) * noise;
B = unit_circle * inner_scale + randn(n,2) * noise;

scatter(A(:,1),A(:,2));
hold on
scatter(B(:,1),B(:,2));

% compute Hausdorff distance by brute force
D = pdist2(A,B);
fhd = max(min(D,[],1));
rhd = max(min(D,[],2));
hd = max(fhd,rhd);

% compute modified Hausdorff distance by brute force
fhd = mean(min(D,[],1));
rhd = mean(min(D,[],2));
mhd = max(fhd, rhd);

% use the mex function
% mhdm = ModHausdorffDistMex(A,B);

% assert(mhd == mhdm);

% now use KDTree
a_kdt = KDTreeSearcher(A);
b_kdt = KDTreeSearcher(B);

[~, fhd] = knnsearch(a_kdt,B,'K',1);
[~, rhd] = knnsearch(b_kdt,A,'K',1);
knn_hd = max(max(fhd), max(rhd));
knn_mhd = max(mean(fhd), mean(rhd));

assert(knn_hd == hd);
assert(knn_mhd == mhd);