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

[tvijlbrief]

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.

[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;

}

[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.

[tvijlbrief]

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

010A00002040710000010000000102000000160000005EBA490CEA0302418D976E9218991C4185EB51B8D403024185EB513810991C41273108ACBC030241295C8F42F3981C414E621058B20302416F128340EB981C4179E92631A8030241D122DB79E5981C4175931804D102024123DBF97E92981C41666666660902024179E9263147981C411904560EC8040241A69BC4A01A981C4166666666CC0402413BDF4F8D1A981C41105839B4D00402417D3F35DE1A981C414260E5D0D4040241B07268911B981C410E2DB29DD8040241C976BE9F1C981C4117D9CEF7DB040241BA490C021E981C41295C8FC2DE040241273108AC1F981C41FED478E9E0040241D34D629021981C414E621058E20402410E2DB29D23981C4175931804E3040241295C8FC225981C41448B6CE7E20402412DB29DEF27981C4100000000E2040241F6285C0F2A981C4108AC1C5AE00402418D976E122C981C41E7FBA9F1BA0402410000000050981C415EBA490CEA0302418D976E9218991C41

010A000020407100000100000001020000001F000000EC51B81EA60A024108AC1CDABD971C41A8C64B37D10A024179E92631BB971C41273108ACD40A0241666666E6BC971C4154E3A59BD60A0241713D0A57C2971C4179E92631D60A024179E926B1C6971C41F4FDD478910A0241000000000F981C41CFF753E3BF0A024177BE9F1A1B981C41E9263108BD0A0241D7A370BD1E981C413BDF4F8D710A0241C520B07267981C414E6210586C0A0241A8C64BB76D981C410C022B876A0A024177BE9F9A70981C411F85EB516A0A0241105839B473981C419EEFA7C6730A0241C74B378984981C41022B8716940A0241C1CAA1C5B8981C41A245B6F3A00A0241E7FBA9F1C4981C4193180456A70A02416891ED7CC8981C41C74B3789AE0A0241560E2D32CB981C417B14AE47EB0A02414260E5D0DC981C41E9263108F10A0241643BDFCFDE981C41F0A7C64BF20A024146B6F3FDDF981C410AD7A370F20A024191ED7C3FE1981C41DD240681F10A024121B072E8E2981C4146B6F3FDEE0A0241B81E856BE5981C41D7A3703DEC0A0241986E1203E8981C41355EBA49E30A02413BDF4F8DF0981C41AAF1D24D890A0241295C8FC2D9981C41D34D6210400A0241448B6C6771981C41F4FDD4788F0A0241931804D622981C41B0726891610A0241FCA9F15216981C41B29DEFA7A80A02411283C0CACB971C41EC51B81EA60A024108AC1CDABD971C41

[pramsey]

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

[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...

[tvijlbrief]

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

[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.

[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.

[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.

[tvijlbrief]

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

[pramsey]

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