Opened 3 years ago

Closed 2 years ago

#4336 closed defect (worksforme)

ST_LineToCurve generates wrong Arcs for Epsilon 1E-1 (+ potential fix)

Reported by: tvijlbrief Owned by: pramsey
Priority: medium Milestone: PostGIS 2.5.3
Component: postgis Version: 2.4.x
Keywords: Cc:

Description

See also https://github.com/tomtor/postgis/blob/master/patch

The two attached images show the issues.

The first image shows that the large light blue arc which clearly exceeds the 0.1m tolerance. This is caused by line 23 in the diff file: comparing an angle difference in radians to the EPSILON value 0.1 (m).

I opted for excepting angles which differ by at most 25% but it would probably be more consistent to compare the distances between a2/a3 and a3/candidate so that they are within (2*) Epsilon.

The second image shows another issue. The current code requires that points which EXTENT the candidate arc have similar angles. However, we do not have such a constraint for a1/a2/a3 which form the ORIGINAL candidate arc. This is inconsistent.

I added a test which require that difference in distance between a1/a2 and a2/a3 is at most 2 * Epsilon.

These errors will be much rarer when using the default 1E-8 epsilon, but the issues are the same.

Attachments (2)

BadCurve.PNG (90.1 KB) - added by tvijlbrief 3 years ago.
BadCurve2.PNG (28.7 KB) - added by tvijlbrief 3 years ago.

Download all attachments as: .zip

Change History (15)

Changed 3 years ago by tvijlbrief

Attachment: BadCurve.PNG added

Changed 3 years ago by tvijlbrief

Attachment: BadCurve2.PNG added

comment:1 Changed 3 years ago by tvijlbrief

I think this is a better alternative for the original angle test (use epsilon):

                int a2_side = lw_segment_side(t1, t3, t2);
                int b_side  = lw_segment_side(t1, t3, tb);
#if 0 // old code
                double angle1 = lw_arc_angle(t1, t2, t3);
                double angle2 = lw_arc_angle(t2, t3, tb);

                /* Is the angle similar to the previous one ? */
                diff = fabs(angle1 - angle2);
                LWDEBUGF(4, " angle1: %g, angle2: %g, diff:%g", angle1, angle2, diff);
                if ( diff > angle1 / 4 || diff > angle2 / 4)
#else // use Epsilon distance
                double distance3 = distance2d_pt_pt(t3, tb);
                diff = fabs(distance3 - distance2);
                if ( diff > 2 * EPSILON_SQLMM )
#endif
                {
                        return LW_FALSE;
                }
Last edited 3 years ago by tvijlbrief (previous) (diff)

comment:2 Changed 3 years ago by robe

Milestone: PostGIS 2.5.2PostGIS 2.5.3

comment:3 Changed 3 years ago by pramsey

I feel like your fix is not necessarily wrong, but neither is it necessarily correct. I see at least two other failure modes in it. I'd be much helped by having the actual input data that cause the incorrect output arcs, so I can check my own approaches to the fix rather than blindly applying something.

Last edited 3 years ago by pramsey (previous) (diff)

comment:4 Changed 2 years ago by tvijlbrief

These are the 2 geometries ( geometry(CurvePolygon?,28992) ):

010A00002040710000010000000102000000160000005EBA490CEA0302418D976E9218991C4185EB51B8D403024185EB513810991C41273108ACBC030241295C8F42F3981C414E621058B20302416F128340EB981C4179E92631A8030241D122DB79E5981C4175931804D102024123DBF97E92981C41666666660902024179E9263147981C411904560EC8040241A69BC4A01A981C4166666666CC0402413BDF4F8D1A981C41105839B4D00402417D3F35DE1A981C414260E5D0D4040241B07268911B981C410E2DB29DD8040241C976BE9F1C981C4117D9CEF7DB040241BA490C021E981C41295C8FC2DE040241273108AC1F981C41FED478E9E0040241D34D629021981C414E621058E20402410E2DB29D23981C4175931804E3040241295C8FC225981C41448B6CE7E20402412DB29DEF27981C4100000000E2040241F6285C0F2A981C4108AC1C5AE00402418D976E122C981C41E7FBA9F1BA0402410000000050981C415EBA490CEA0302418D976E9218991C41

010A000020407100000100000001020000001F000000EC51B81EA60A024108AC1CDABD971C41A8C64B37D10A024179E92631BB971C41273108ACD40A0241666666E6BC971C4154E3A59BD60A0241713D0A57C2971C4179E92631D60A024179E926B1C6971C41F4FDD478910A0241000000000F981C41CFF753E3BF0A024177BE9F1A1B981C41E9263108BD0A0241D7A370BD1E981C413BDF4F8D710A0241C520B07267981C414E6210586C0A0241A8C64BB76D981C410C022B876A0A024177BE9F9A70981C411F85EB516A0A0241105839B473981C419EEFA7C6730A0241C74B378984981C41022B8716940A0241C1CAA1C5B8981C41A245B6F3A00A0241E7FBA9F1C4981C4193180456A70A02416891ED7CC8981C41C74B3789AE0A0241560E2D32CB981C417B14AE47EB0A02414260E5D0DC981C41E9263108F10A0241643BDFCFDE981C41F0A7C64BF20A024146B6F3FDDF981C410AD7A370F20A024191ED7C3FE1981C41DD240681F10A024121B072E8E2981C4146B6F3FDEE0A0241B81E856BE5981C41D7A3703DEC0A0241986E1203E8981C41355EBA49E30A02413BDF4F8DF0981C41AAF1D24D890A0241295C8FC2D9981C41D34D6210400A0241448B6C6771981C41F4FDD4788F0A0241931804D622981C41B0726891610A0241FCA9F15216981C41B29DEFA7A80A02411283C0CACB971C41EC51B81EA60A024108AC1CDABD971C41

Last edited 2 years ago by tvijlbrief (previous) (diff)

comment:5 Changed 2 years ago by pramsey

Sorry, and with these geometries, you are just running a st_linetocurve(st_curvetoline()) cycle on them?

comment:6 Changed 2 years ago by tvijlbrief

I am not sure what you mean. These two geometries are part of a larger collection of CurvePolygons? (the black lines in the two images) which I select into a new relation with select(st_linetocurve(org_geom)).

The two images show the visualisation (in light blue) of the result in QGIS. I do not use a st_curvetoline() myself, but QGIS might do this when drawing the Postgis geometrie on the canvas...

Last edited 2 years ago by tvijlbrief (previous) (diff)

comment:7 Changed 2 years ago by tvijlbrief

Note that I only show the segments which are CURVES in the light blue result set.

comment:8 Changed 2 years ago by pramsey

OK, now I'm quite confused. The geometries you are giving are, in WKT, this

CURVEPOLYGON((147581.256 468550.143,147578.59 468548.055,147575.584 468540.815,147574.293 468538.813,147573.024 468537.369,147546.127 468516.624,147521.175 468497.798,147609.007 468486.657,147609.55 468486.638,147610.088 468486.717,147610.602 468486.892,147611.077 468487.156,147611.496 468487.502,147611.845 468487.918,147612.114 468488.391,147612.293 468488.904,147612.377 468489.44,147612.363 468489.984,147612.25 468490.515,147612.044 468491.018,147607.368 468500,147581.256 468550.143))

and

CURVEPOLYGON((147796.765 468463.463,147802.152 468462.798,147802.584 468463.225,147802.826 468464.585,147802.774 468465.673,147794.184 468483.75,147799.986 468486.776,147799.629 468487.685,147790.194 468505.862,147789.543 468507.429,147789.316 468508.151,147789.29 468508.926,147790.472 468513.134,147794.511 468526.193,147796.119 468529.236,147796.917 468530.122,147797.817 468530.799,147805.41 468535.204,147806.129 468535.703,147806.287 468535.998,147806.305 468536.312,147806.188 468536.727,147805.874 468537.355,147805.53 468538.003,147804.411 468540.138,147793.163 468534.44,147784.008 468508.351,147793.934 468488.709,147788.196 468485.581,147797.082 468466.948,147796.765 468463.463))

But these geometries don't actually include curves... the WKT for CURVEPOLYGON requires that the curved portions be appropriately tagged as CIRCULARSTRING. Untagged portions are considered to be ordinary LINEARRINGs. So now I'm really not sure what's the problem precisely. The picture isn't what I want.

I'd like to see:

  • An input LINESTRING.
  • The simple, three-point CIRCULARSTRING you'd expect to get from ST_LineToCurve .
  • The CIRCULARSTRING you *do* get from ST_LineToCurve

Basically, I need to be able to replicate your condition, on my own machine. That way I can (a) verify that your fix works for me and (b) verify if I have a different fix that also works.

Last edited 2 years ago by pramsey (previous) (diff)

comment:9 Changed 2 years ago by tvijlbrief

Yes, they do not include curves. You can change the signature to linestring, that is not relevant.

Just feed the first geometry to st_linetocurve() from an unchanged postgis with epsilon 1e-1 and check the result. It will generate a result linestring with an error epsilon > 10m.

comment:10 Changed 2 years ago by pramsey

ST_LineToCurve() generates curves from lines that are stroked from curves (collections of linear segments that can be approximated with a similar curve. When I run:

select st_astext(st_linetocurve('LINESTRING(147581.256 468550.143,147578.59 468548.055,147575.584 468540.815,147574.293 468538.813,147573.024 468537.369,147546.127 468516.624,147521.175 468497.798,147609.007 468486.657,147609.55 468486.638,147610.088 468486.717,147610.602 468486.892,147611.077 468487.156,147611.496 468487.502,147611.845 468487.918,147612.114 468488.391,147612.293 468488.904,147612.377 468489.44,147612.363 468489.984,147612.25 468490.515,147612.044 468491.018,147607.368 468500,147581.256 468550.143)'::geometry))

I get back

LINESTRING(147581.256 468550.143,147578.59 468548.055,147575.584 468540.815,147574.293 468538.813,147573.024 468537.369,147546.127 468516.624,147521.175 468497.798,147609.007 468486.657,147609.55 468486.638,147610.088 468486.717,147610.602 468486.892,147611.077 468487.156,147611.496 468487.502,147611.845 468487.918,147612.114 468488.391,147612.293 468488.904,147612.377 468489.44,147612.363 468489.984,147612.25 468490.515,147612.044 468491.018,147607.368 468500,147581.256 468550.143)

which is the same as the input. The function is saying there's not enough similarity in the edges to approximate a curve so it's returning the input untouched. The ST_LineToCurve() also doesn't take an epsilon, though ST_CurveToLine() does.

Can you summarize your issue with the smallest possible collection of SQL and data? I can't keep tossing darts at the board and hoping I get the right combination.

comment:11 Changed 2 years ago by tvijlbrief

You have to change the Postgis COMPILE epsilon to 1e-1 to reproduce the error because st_linetocurve() does not have that parameter.

comment:12 Changed 2 years ago by pramsey

Again, I'm not seeing any difference, even changing the compile value.

comment:13 Changed 2 years ago by pramsey

Resolution: worksforme
Status: newclosed
Note: See TracTickets for help on using tickets.