5

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)

  • What version of GEOS are you using, paste the output of SELECT postgis_full_version(); – Evan Carroll Aug 23 '17 at 19:14
  • @EvanCarroll Please see the update for the PostGIS/GEOS versions. – tinlyx Aug 23 '17 at 19:19
  • I'm confused at the problem now. You're using PostGIS 2.3.3 why don't you just use the hausdorffdistance function in PostGIS? I don't get what you're trying to do. – Evan Carroll Aug 23 '17 at 19:21
  • @EvanCarroll The Hausdorff distance is a real distance, and symmetric. I only needed the one-sided Hausdorff distance, which is asymmetric. Please see e.g.: cgm.cs.mcgill.ca/~godfried/teaching/cg-projects/98/normand/… for definitions – tinlyx Aug 23 '17 at 19:23
2

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,

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.

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.

  • Thanks for the pointer. In DiscreteHausdorffDistance.h: double orientedDistance() { computeOrientedDistance(g0, g1, ptDist); return ptDist.getDistance(); } seem interesting. Do you know what it is and how to expose it to PostGIS for testing? – tinlyx Aug 23 '17 at 19:48
  • @tinlyx updated. – Evan Carroll Aug 23 '17 at 22:56
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)

  • Thanks. Could you elaborate the process a bit or provide a link? I'm not very familiar with it. – tinlyx Aug 23 '17 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.