Intersection tests with curved polygons

Wed Jan 8 01:13:56 PST 2025

Hi Paul,
apologies for the confusion, I just can't type the right word sometimes.
Only ST_Intersects is being used.
Here are again the numbers with the actual query (reminder, it's being run
against the original table, hence the different table name compared to
the test I shared):

Baseline (intersection), provides the wrong answer (two polygons):

SELECT ogc_fid FROM kiinteisto_alue WHERE geom && ST_GeomFromText('POINT
(25492818 6677399.98)', 3879) and ST_Intersects(geom,
ST_GeomFromText('POINT (25492818 6677399.98)', 3879));
tps = 10227.810087 ((without initial connection time)

Distance, provides the right answer:

SELECT ogc_fid FROM kiinteisto_alue WHERE geom && ST_GeomFromText('POINT
(25492818 6677399.98)', 3879) and ST_Distance(geom, ST_GeomFromText('POINT
(25492818 6677399.98)', 3879)) = 0;
tps = 19989.627257 (without initial connection time)

Manual linearization with higher accuracy (provides the right answer):

SELECT ogc_fid FROM kiinteisto_alue WHERE geom && ST_GeomFromText('POINT
(25492818 6677399.98)', 3879) and ST_Intersects(ST_CurveToLine(geom, 0.01,
1, 1), ST_GeomFromText('POINT (25492818 6677399.98)', 3879))
tps = 8984.594159 (without initial connection time)

In the full dataset, there are around 10k records having actual curved
edges, out of 36k.
I was guessing the ideal would be to check if curves are there, and change
on the fly the test.
Something like this:

SELECT ogc_fid FROM kiinteisto_alue WHERE geom && ST_GeomFromText('POINT
(25492818 6677399.98)', 3879) and
case
  when ST_HasArc(geom) then ST_Distance(geom, ST_GeomFromText('POINT
(25492818 6677399.98)', 3879)) = 0
  else ST_Intersects(geom, ST_GeomFromText('POINT (25492818 6677399.98)',
3879))
end;

And yet, a pgbench of the above only returns 14500 TPS, whilst the simple
distance test for all,
gets up to around 20000.... looks like a ST_Distance test is faster even
for geometries with straight segments.... isn't this odd?
In my mind, checking if two geometries interfere (in any way) should be
less expensive than calculating their
exact distance.

By the way, thanks for the details, I love to hear about how things are
implemented!

Regards,

Andrea Aime

==

GeoServer Professional Services from the experts!

Visit http://bit.ly/gs-services-us for more information.

==

Ing. Andrea Aime
@geowolf
Technical Lead

GeoSolutions Group
phone: +39 0584 962313

fax:     +39 0584 1660272

mob:   +39  339 8844549

https://www.geosolutionsgroup.com/

http://twitter.com/geosolutions_it

-------------------------------------------------------

Con riferimento alla normativa sul trattamento dei dati personali (Reg. UE
2016/679 - Regolamento generale sulla protezione dei dati “GDPR”), si
precisa che ogni circostanza inerente alla presente email (il suo
contenuto, gli eventuali allegati, etc.) è un dato la cui conoscenza è
riservata al/i solo/i destinatario/i indicati dallo scrivente. Se il
messaggio Le è giunto per errore, è tenuta/o a cancellarlo, ogni altra
operazione è illecita. Le sarei comunque grato se potesse darmene notizia.

This email is intended only for the person or entity to which it is
addressed and may contain information that is privileged, confidential or
otherwise protected from disclosure. We remind that - as provided by
European Regulation 2016/679 “GDPR” - copying, dissemination or use of this
e-mail or the information herein by anyone other than the intended
recipient is prohibited. If you have received this email by mistake, please
notify us immediately by telephone or e-mail

On Tue, Jan 7, 2025 at 7:36 PM Paul Ramsey <pramsey at cleverelephant.ca>
wrote:

> I’m not 100% sure this is apples and apples?
>
> On Jan 7, 2025, at 2:55 AM, Andrea Aime <andrea.aime at geosolutionsgroup.com>
> wrote:
>  testing 3 scenarios:
>
>    - The current query using ST_Intersection as is (returns 2 polygons,
>    wrong result but I'm treating it as the baseline)
>
> ST_Intersection() is an expensive operation, constructing new geometry.
> ST_Intersects() is a cheaper operation, returning true/false.
>
>
>    - Using && + ST_Distance
>
> ST_Distance() is cheaper than ST_Intersection(), and similar to
> ST_Intersects() in how it calculates an answer.
>
>
>    - Using && + linearized intersection, with a precision good enough to
>    get the correct result (0.01 meters)
>
> Still using Intersection() so on an expensive path.
>
> Some things to note in general:
>
> - ST_Intersection() will always delegate to GEOS and if handed curves will
> linearize them first. Delegating to GEOS does have some fixed overheads, a
> full copy of the geometry is made, and GEOS set up.
> - ST_Intersects() will sometimes delegate to GEOS, one of those cases
> being handed a curve polygon. Which will again coast you linearization and
> a full copy.
> - Finer linearization will naturally result in more segments and more
> processing cost.
> - ST_Distance() with native curve support trades the expense of some extra
> trig doing edge/arc calculations for avoiding the trip to GEOS and
> processing a much smaller number of edges than the linearized path.
>
> I find the results surprsing:
>
> > pgbench -U cite helsinki -f baseline.sql  -c 32 -t 2000
> tps = 10227.810087 ((without initial connection time)
>
> > pgbench -U cite helsinki -f distance.sql  -c 32 -t 2000
> tps = 19989.627257 (without initial connection time)
>
> > pgbench -U cite helsinki -f linearized.sql  -c 32 -t 2000
> tps = 8984.594159 (without initial connection time)
>
> So... testing with the distance is twice as fast as the other options?
> Wow, have we been doing intersection tests wrong all this time? ROFL
> Other possible ideas:
>
>    - There is something specific to having curves in the mix, and having
>    to pay the cost of linearization makes distance competitive?
>    - The specific dataset is playing an important role in the result
>
> I am not, in general, surprised at your results, knowing the different
> code paths your tests could be running through.
>
> P.
>
>
> Ah, since there seems to be a real issue, I've also opened a ticket in
> trac: https://trac.osgeo.org/postgis/ticket/5832#ticket
>
> Regards,
> Andrea Aime
>
> ==
>
> GeoServer Professional Services from the experts!
> Visit http://bit.ly/gs-services-us for more information.
>
> ==
>
> Ing. Andrea Aime
> @geowolf
> Technical Lead
>
> GeoSolutions Group
> phone: +39 0584 962313
> fax:     +39 0584 1660272
> mob:   +39  339 8844549
>
> https://www.geosolutionsgroup.com/
> http://twitter.com/geosolutions_it
> -------------------------------------------------------
>
> Con riferimento alla normativa sul trattamento dei dati personali (Reg. UE
> 2016/679 - Regolamento generale sulla protezione dei dati “GDPR”), si
> precisa che ogni circostanza inerente alla presente email (il suo
> contenuto, gli eventuali allegati, etc.) è un dato la cui conoscenza è
> riservata al/i solo/i destinatario/i indicati dallo scrivente. Se il
> messaggio Le è giunto per errore, è tenuta/o a cancellarlo, ogni altra
> operazione è illecita. Le sarei comunque grato se potesse darmene notizia.
>
> This email is intended only for the person or entity to which it is
> addressed and may contain information that is privileged, confidential or
> otherwise protected from disclosure. We remind that - as provided by
> European Regulation 2016/679 “GDPR” - copying, dissemination or use of this
> e-mail or the information herein by anyone other than the intended
> recipient is prohibited. If you have received this email by mistake, please
> notify us immediately by telephone or e-mail
>
>
> On Sun, Dec 22, 2024 at 4:44 PM Andrea Aime <
> andrea.aime at geosolutionsgroup.com> wrote:
>
>> Hi Paul,
>> thanks a lot for following up. Comments inline below.
>>
>> These are literally CurvePolygon type?
>>>
>>
>> The column type is just "geometry(Geometry,3879)", while ST_GeometryType
>> returns "multisurface" for both.
>> When doing a ST_AsText instead, you'll get something like:
>>
>> MULTISURFACE(CURVEPOLYGON(COMPOUNDCURVE((...
>>
>> for both.
>>
>>
>>> It’s probably getting caught in our lack of full curve support.
>>> I would be interested in the ST_Distance between the point and those two
>>> CurvePolygons. (Because, for distance, we have a postgis-native
>>> implementation that supports curves).
>>>
>>
>> =# SELECT ogc_fid, ST_Distance(ST_GeomFromText('POINT (25492818
>> 6677399.98)', 3879), geom) FROM testdata;
>>
>>  ogc_fid |     st_distance
>> ---------+---------------------
>>     1258 | 0.01234572446598792
>>    12875 |                   0
>> (2 rows)
>>
>> Indeed, the correct answer, 12875 contains the point, while the other
>> polygon is close to it.
>>
>>
>>> Whereas for intersection, the calculation is delegated to GEOS *after
>>> linearizing the inputs*. In that linearization, could sit the logically
>>> problem you’re seeing.
>>>
>>
>> Let's check with different tolerances... yes, changing the tolerance
>> changes the result:
>>
>> =#  SELECT ogc_fid FROM testdata WHERE ST_Intersects(ST_CurveToLine(geom,
>> 0.01, 1, 1), ST_GeomFromText('POINT (25492818 6677399.98)', 3879));
>>  ogc_fid
>> ---------
>>    12875
>> (1 row)
>>
>> #  SELECT ogc_fid FROM testdata WHERE ST_Intersects(ST_CurveToLine(geom,
>> 0.02, 1, 1), ST_GeomFromText('POINT (25492818 6677399.98)', 3879));
>>  ogc_fid
>> ---------
>>     1258
>> (1 row)
>>
>> In the immediate future, I guess I could have the GeoTools PostGIS store
>> use either approach, when knowing curves are involved...
>> First using && to perform a first rough filter, and then either use either
>> * ST_Distance equals to 0
>> * An explicit linearization with a target tolerance (this is an urban
>> application, so I'm guessing they will need centimeter, if not millimeter,
>> precision)
>> .
>> Is there a clear winner here in terms of performance, or performance of
>> distance vs linearized intersection is more contextual to the geometries
>> involved?
>>
>> Cheers
>> Andrea
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250108/863b3606/attachment.htm>