Directed Hausdorff distance in PostGIS?

Question

I was wondering how to implement a directed (or one-sided) Hausdorff distance function in PostGIS?

I looked into the PostGIS 2.3.3 source code, and it seems that ST_HausdorffDistance is implemented as a simple C wrapper function that calls the corresponding GEOS function as follows:

PG_FUNCTION_INFO_V1(hausdorffdistance);
Datum hausdorffdistance(PG_FUNCTION_ARGS)
{
    GSERIALIZED *geom1;
    GSERIALIZED *geom2;
    GEOSGeometry *g1;
    GEOSGeometry *g2;
    double result;
    int retcode;

    POSTGIS_DEBUG(2, "hausdorff_distance called");

    geom1 = PG_GETARG_GSERIALIZED_P(0);
    geom2 = PG_GETARG_GSERIALIZED_P(1);

    if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
        PG_RETURN_NULL();

    initGEOS(lwpgnotice, lwgeom_geos_error);

    g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
    if ( 0 == g1 )   /* exception thrown at construction */
    {
        HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
        PG_RETURN_NULL();
    }

    g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
    if ( 0 == g2 )   /* exception thrown */
    {
        HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
        GEOSGeom_destroy(g1);
        PG_RETURN_NULL();
    }

    retcode = GEOSHausdorffDistance(g1, g2, &result);
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);

    if (retcode == 0)
    {
        HANDLE_GEOS_ERROR("GEOSHausdorffDistance");
        PG_RETURN_NULL(); /*never get here */
    }

    PG_FREE_IF_COPY(geom1, 0);
    PG_FREE_IF_COPY(geom2, 1);

    PG_RETURN_FLOAT8(result);
}

I don't have much experience with GEOS or interfacing it with PostGIS. The only thing I found in GEOS is that in that in geos::algorithm::distance::DiscreteHausdorffDistance, it is said that there is an double orientedDistance () member function. I am not sure if the documentation is even current.

Can some with the knowledge suggest how to adapt the above code to obtain the directed (or one-sided) Hausdorff distance from GEOS?

--EDIT--

The PostGIS info from SELECT postgis_full_version();' is as below (Ubuntu 16.04).

POSTGIS="2.3.3 r15473" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" (core procs from "2.3.2 r15302" need upgrade) TOPOLOGY (topology procs from "2.3.2 r15302" need upgrade) RASTER (raster procs from "2.3.2 r15302" need upgrade) (sfcgal procs from "2.3.2 r15302" need upgrade)

Answer 1

I'm not sure what you're talking about but the GEOS implimentation is pretty straight forward,

int GEOS_DLL GEOSHausdorffDistance(
   const GEOSGeometry *g1,
   const GEOSGeometry *g2,
   double *dist
);  

That is to say, it takes geometry's and mutates a pointer to a double returning an int that represents status. There are no extra arguments. There are no methods of telling it "asymmetric" or "symmetric" or "one-sided". So fullstop, if you're using GEOS you're using, DiscreteHausdorffDistance.cpp

You can see the implementation here,

• DiscreteHausdorffDistance.h
• DiscreteHausdorffDistance.cpp

Update

I think you're way over your head with this, but here is the breakdown of the source.

PostGIS Bindings

There are two functions in PostGIS, you can see them made available to the user here

1. ST_HausdorffDistance(geom1 geometry, geom2 geometry)

   1. which dispatches to hausdorffdistance.
   2. hausdorffdistance references GEOSHausdorffDistance.

2. ST_HausdorffDistance(geom1 geometry, geom2 geometry, float8)

   1. which dispatches to hausdorffdistancedensify.
   2. hausdorffdistancedensify function references GEOSHausdorffDistanceDensify.

LibGEOS

LibGEOS provides two corresponding functions

1. GEOSHausdorffDistance(const Geometry *g1, const Geometry *g2, double *dist)

   1. GEOSHausdorffDistance dispatches to GEOSHausdorffDistance_r
   2. GEOSHausdorffDistance_r's implementation ultimately invokes DiscreteHausdorffDistance::distance(*g1, *g2)

2. GEOSHausdorffDistanceDensify(const Geometry *g1, const Geometry *g2, double densifyFrac, double *dist)

   1. GEOSHausdorffDistanceDensify dispatches to GEOSHausdorffDistanceDensify_r
   2. GEOSHausdorffDistanceDensify_r's implementation ultimately invokes DiscreteHausdorffDistance::distance(*g1, *g2, densifyFrac)

DiscreteHausdorffDistance::distance

Now you can see all roads end up at DiscreteHausdorffDistance::distance which is overloaded.

• 2 argument: double DiscreteHausdorffDistance::distance(const geom::Geometry& g0, const geom::Geometry& g1)
• 3 argument: double DiscreteHausdorffDistance::distance(const geom::Geometry& g0, const geom::Geometry& g1, double densifyFrac)

Both of those

1. create a DiscreteHausdorffDistance
2. Call distance on the object they create.
3. distance calls computeOrientedDistance(g0, g1, ptDist) and computeOrientedDistance(g1, g0, ptDist) which is defined here
4. computeOrientedDistance mutates a private ptDist by calling ptDist.setMaximum

So the answer to your question: computeOrientedDistance is already called by both ST_HausdorffDistance. If you want a directional version, like computeOrientedDistance I think you'll have to pull the above apart and patch it in LibGEOS via a new DiscreteHausdorffDistance::orientedDistance and then into PostGIS. Or draw up a test suite and pay a consultant. While there are probably a few other software programmers here that could do it, it's way out of the expertise of this community (Paul Ramsey being an obvious exception, etc).

Generally, we don't implement software suggestions on GIS.SE. At least, I don't see that activity if it happens. We answer them with software already developed. The short answer here is that it's not done.

Answer 2

You need to first expose the "oriented" as a new CAPI function, then you would use the new function from postgis (will depend on new GEOS release)