Hausdorff Distance values odd

[plooney]

I'm using a python wrapper I wrote to test some functions. I expect the Hausdorff Distance to be close to 0 but I get nonzero values in Rust.Geo. Converting the shapes to shapely gives me a value close to zero.

If I simplify the multi polygon I can get very different values.

!uv pip install -U rust_geo_python
from rust_geo_python import RustPolygon, union
from shapely import from_wkt

def make_circle(centre, radius, nVerts=200):
    direction = -1

    pts = []
    for v in range(nVerts):
        theta = direction*v*2*np.pi/nVerts
        pts.append([radius*np.cos(theta), radius*np.sin(theta)])
    pts = pts + np.asarray(centre)
    return np.asarray(pts)

def test_buffering():
    rect_verts_ext = np.array([[2,2],[2,-2],[-2,-2],[-2,2]]).astype(float)
    rect_verts_int = np.array([[1,1],[1,-1],[-1,-1],[-1,1]]).astype(float)[::-1]
    rect_polygon = RustPolygon(rect_verts_ext, [rect_verts_int])
    rect_verts_int = np.array([[0.5,0.5],[0.5,-0.5],[-0.5,-0.5],[-0.5,0.5]]).astype(float)[::-1]
    rect_polygon_small_int = RustPolygon(rect_verts_ext, [rect_verts_int])
    n_points = 100
    x = np.linspace(-2,2,n_points)
    y = np.linspace(-2,2,n_points)
    left_points = list(zip(np.array([-2]*n_points), y))
    right_points = list(zip(np.array([2]*n_points), y))
    bottom_points = list(zip(x, [-2]*n_points))
    top_points = list(zip(x, [2]*n_points))
    circles = [RustPolygon(make_circle(x, 0.5, nVerts=100),[]) for x in left_points + right_points + top_points + bottom_points]
    rect_union = union(circles + [rect_polygon, rect_polygon_small_int])
    buffered_rect = rect_polygon.buffer(0.5)
    print(f"Rust value {rect_union.hausdorff_distance(buffered_rect)}")
    return buffered_rect, rect_union, rect_polygon_small_int

rect_buffer, rect_union, rect_small = test_buffering()

import shapely
print(f"Shapely value {shapely.hausdorff_distance(from_wkt(rect_buffer.to_wkt()), from_wkt(rect_union.to_wkt()))}")

[kylebarron]

I'm not familiar with this area of the geo code myself, but I think you're most likely to get a helpful response if you add a reproducible example in rust, not Python

[plooney]

@kylebarron as suggested.

use geo::{HausdorffDistance, Simplify, Validation};
use std::f64::consts::PI;

use geo::{AffineOps, AffineTransform, unary_union};
use geo::{BooleanOps, Buffer, LineString, MultiPolygon, Point, Polygon, Rect, coord};
use ndarray::linspace;
use wkt::ToWkt;

fn main() {
    println!("Hello, world!");

    let buffered_rect = make_rect_poly_buffered();
    let rect_small_int = make_rect_poly();
    let circles = make_circles(200);
    let union = unary_union(&circles).union(&rect_small_int);
    let hd = union.hausdorff_distance(&buffered_rect);
    println!("Expect close to 0. Got {}", hd);
    println!();
    println!("{}", union.is_valid());
    println!();
    println!("{}", buffered_rect.wkt_string());

    //Oddly simplifying gets an answer closer to zero

    let simple_union = union.simplify(0.01);
    let hd = simple_union.hausdorff_distance(&buffered_rect);
    println!("Expect close to 0. Got {}", hd);
    println!();
    println!("{}", simple_union.wkt_string());
    println!();
    println!("{}", buffered_rect.wkt_string());
}

fn make_rect_poly() -> Polygon<f64> {
    let (x1, y1) = (-2.0, -2.0) as (f64, f64);
    let (x2, y2) = (2.0, 2.0) as (f64, f64);
    let rect_ext = Rect::new(coord! {x: x1, y: y1}, coord! {x: x2, y: y2 });
    let (x1, y1) = (-0.5, -0.5) as (f64, f64);
    let (x2, y2) = (0.5, 0.5) as (f64, f64);
    let rect_int = Rect::new(coord! {x: x1, y: y1}, coord! {x: x2, y: y2 });
    Polygon::new(
        rect_ext.to_polygon().exterior().clone(),
        vec![rect_int.to_polygon().exterior().clone()],
    )
}

fn make_rect_poly_buffered() -> MultiPolygon<f64> {
    let (x1, y1) = (-2.0, -2.0) as (f64, f64);
    let (x2, y2) = (2.0, 2.0) as (f64, f64);
    let rect_ext = Rect::new(coord! {x: x1, y: y1}, coord! {x: x2, y: y2 });
    let (x1, y1) = (-1.0, -1.0) as (f64, f64);
    let (x2, y2) = (1.0, 1.0) as (f64, f64);
    let rect_int = Rect::new(coord! {x: x1, y: y1}, coord! {x: x2, y: y2 });
    Polygon::new(
        rect_ext.to_polygon().exterior().clone(),
        vec![rect_int.to_polygon().exterior().clone()],
    )
    .buffer(0.5)
}

fn make_circles(n_points: usize) -> Vec<Polygon<f64>> {
    let (x1, y1) = (-2.0, -2.0) as (f64, f64);
    let (x2, y2) = (2.0, 2.0) as (f64, f64);
    let xs = linspace(x1, x2, n_points);
    let ys = linspace(y1, y2, n_points);
    let top_circles = xs
        .map(|x| make_circle(x, 2.0, 0.5, 200))
        .collect::<Vec<Polygon<f64>>>();
    let xs = linspace(x1, x2, n_points);
    let bottom_circles = xs
        .map(|x| make_circle(x, -2.0, 0.5, 200))
        .collect::<Vec<Polygon<f64>>>();
    let left_circles = ys
        .map(|x| make_circle(-2.0, x, 0.5, 200))
        .collect::<Vec<Polygon<f64>>>();
    let ys = linspace(x1, x2, n_points);
    let right_circles = ys
        .map(|x| make_circle(2.0, x, 0.5, 200))
        .collect::<Vec<Polygon<f64>>>();
    [top_circles, bottom_circles, left_circles, right_circles].concat()
}

fn make_circle(x: f64, y: f64, radius: f64, n: i32) -> Polygon<f64> {
    let direction = -1.0;
    let points = (0..n)
        .map(|i| {
            let i = i as f64;
            let n = n as f64;
            let theta = direction * i * 2.0 * PI / n;
            let x = radius * f64::cos(theta);
            let y = radius * f64::sin(theta);
            Point::new(x, y)
        })
        .collect::<LineString<f64>>();
    let polygon = Polygon::new(points, vec![]);
    let transform = AffineTransform::translate(x, y);
    polygon.affine_transform(&transform)
}

[michaelkirk]

You're talking about Hausdorff distance, but your example code also includes buffering. At a glance, I can't tell if this is just another case of #1455

Try rendering the intermediate steps to see if you can answer your own question here.

If you still think there's an issue with the hausdorff implementation, please give your example inputs and the value you expect.

[plooney]

Hi Michael, I've written it out to wkt read it in with shapely and the answer for the non simplified polygon is closer to zero. 2 is not correct.

…

On Mon, 1 Dec 2025, 20:01 Michael Kirk, ***@***.***> wrote: *michaelkirk* left a comment (georust/geo#1456) <#1456 (comment)> You're talking about Hausdorff distance, but your example code also includes buffering. At a glance, I can't tell if this is just another case of #1455 <#1455> Try rendering the intermediate steps to see if you can answer your own question here. If you still think there's an issue with the hausdorff implementation, please give your example inputs and the value you expect. — Reply to this email directly, view it on GitHub <#1456 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABKFSASLWNX7FRVPBSWLIVT37SNBZAVCNFSM6AAAAACNLYLMUWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTKOJYGYYDKNZZGA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[michaelkirk]

I'm not in denial that there is a bug, but could you please provide a simpler reproduction so that we volunteers can focus on fixing the bug rather than trying to understand your example?

[plooney]

I was writing a test using the Hausdorff distance to check that buffering was working as I expected. The two objects should be the same. I suspect that the issue is due to the number of points on the boundary so currently this is the only way I am aware of for reproducing this. I'll have a look in a few days if I can reproduce this in less code.

[urschrei]

Is this is because we're only computing vertex-to-vertex distances, not vertex-to-edge distances, i.e. what we're calculating is discrete Hausdorff distance, and we have in fact documented it as such though the trait method doesn't give a hint: It measures the maximum euclidean distance of a point in one set to the nearest point in another set.

This currently fails:

    #[test]
    fn test_hausdorff_line_edge_distance() {
        use crate::{HausdorffDistance, LineString, coord};

        let a: LineString = vec![coord! {x: 0., y: 0.}, coord! {x: 2., y: 0.}].into();
        let b: LineString = vec![
            coord! {x: 0., y: 0.},
            coord! {x: 1., y: 0.},
            coord! {x: 2., y: 0.},
        ]
        .into();

        let hd = a.hausdorff_distance(&b);
        assert!(hd < 1e-10, "collinear segments: expected HD = 0, got {hd}");
    }

So we need to establish whether the reported bug above is a manifestation of this (I suspect it is). Switching to point-line distance in our algorithm should be easy enough for testing purposes.

[plooney]

Thanks @urschrei that is interesting.

I also suspect that this bug is a manifestation of this. I've not checked but the dense union is going to have some points midway between -2,2 so getting an answer close to 2 makes sense rather than 0 and then getting close to 0 when simplifying when these midpoints are removed.

It seems a little odd that this is the implementation for a Polygon or LineString to me. It should be using the point-line distance in the algorithm.

[michaelkirk]

Yeah, I think it's worth revisiting at least our documentation, if not also our naming conventions. And maybe also making it easier to densify your geometries, or at least documenting that.

Does Point x Line distance make sense for hausdorff distance? I'm not super fluent, but my understanding is that it's an algorithm based on points.

In case you weren't aware @plooney:

1. It sounds like Euclidean.distance would behave more like you're expecting. Is there a reason you're using Hausdorff?

2. For comparing two polygons for similarity, I'll often do something based on xor like: a.xor(b).unsigned_area() / a.unsigned_area()

[urschrei]

On the one hand, our Distance trait should be consistent and I assumed it was always minimum distance no matter the metric space. On the other hand, the semantics are documented.

I think I come down on the side of consistency, even though that'll necessitate a breaking change and leaves open the question of where to integrate discrete distance. Let's see what others have to say.

[urschrei]

Ah, didn't see you'd already commented @michaelkirk, sorry!

[plooney]

Formerly I have worked in medical imaging where Hausdorff distance is a standard measure of similarity of shapes. The Hausdorff distance of a point to a line should be the minimum distance from the point to the line not the minimum distance of the point to one of the vertices.

@michaelkirk

Euclidean distance won't be a good measure of similarity since if there is any overlap it is zero

I ended up doing an area comparison for my single test. What you propose is 1- the intersection/union = 1 - the Jaccard Index I think. That would be fine.

I am not heavily using Hausdorff distance. I just stumbled upon this using it for comparison in a test case.I am just a little troubled by the behaviour though. It seems incorrect. If I wanted to find the discrete Hausdorff distance as is the current behaviour I would convert it to a multi point and then calculate the Hausdorff distance.

[michaelkirk]

As I said, I'm not a big user of Hausdorff distance, so I'm not familiar with its conventions.

Can you point me to an implementation of continuous Hausdorff? e.g. shapely/geos/jts also use discrete Hausdorff distance, though they make it easier to "densify"

I'm doing a bit of research just now. There are some implementations which have a more (IMO) intuitive phrasing of "maximum error tolerance", but I think ultimately it amounts to how much sampling you want to do for what is still a discrete Hausdorff calculation.

I'm definitely in favor of aligning expectations to avoid confusion.

[plooney]

This is a bit more of a can of worms than I anticipated.I understand mathematically how it is defined.

I assumed that the shapely implementation was continuous but I see, as you say, that they densify and then do it discreetly.

It is 3D but VTK sounds like it does it.

https://vtk.org/doc/nightly/html/classvtkHausdorffDistancePointSetFilter.html

"Compute Hausdorff distance between two point sets.

This class computes the relative and hausdorff distances from two point sets (input port 0 and input port 1). If no topology is specified (ie. vtkPointSet or vtkPolyData without vtkPolys), the distances are computed between point location. If polys exist (ie triangulation), the TargetDistanceMethod allows for an interpolation of the cells to ensure a better minimal distance exploration."