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
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,
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
ST_HausdorffDistance(geom1 geometry, geom2 geometry)
- which dispatches to
hausdorffdistance
. hausdorffdistance
referencesGEOSHausdorffDistance
.
- which dispatches to
ST_HausdorffDistance(geom1 geometry, geom2 geometry, float8)
- which dispatches to
hausdorffdistancedensify
. hausdorffdistancedensify
function referencesGEOSHausdorffDistanceDensify
.
- which dispatches to
LibGEOS
LibGEOS provides two corresponding functions
GEOSHausdorffDistance(const Geometry *g1, const Geometry *g2, double *dist)
GEOSHausdorffDistance
dispatches toGEOSHausdorffDistance_r
GEOSHausdorffDistance_r
's implementation ultimately invokesDiscreteHausdorffDistance::distance(*g1, *g2)
GEOSHausdorffDistanceDensify
dispatches toGEOSHausdorffDistanceDensify_r
GEOSHausdorffDistanceDensify_r
's implementation ultimately invokesDiscreteHausdorffDistance::distance(*g1, *g2, densifyFrac)
Now you can see all roads end up at DiscreteHausdorffDistance::distance
which is overloaded.
Both of those
- create a DiscreteHausdorffDistance
- Call
distance
on the object they create. distance
callscomputeOrientedDistance(g0, g1, ptDist)
andcomputeOrientedDistance(g1, g0, ptDist)
which is defined herecomputeOrientedDistance
mutates a privateptDist
by callingptDist.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.
No comments:
Post a Comment