ST_curveToLine bad conversion for some arcs

[claudeunger]

this arc is badly converting : CIRCULARSTRING(2047538.600 7268770.435,2047538.598 7268770.435,2047538.596 7268770.436)

select ST_asText(ST_curveToLine(ST_geomFromText('CIRCULARSTRING(2047538.600 7268770.435,2047538.598 7268770.435,2047538.596 7268770.436)')));

results in : LINESTRING(2047538.6 7268770.435, 2047538.70324625 7268770.45039807, 2047538.80712369 7268770.46071155, 2047538.91138205 7268770.46591557, 2047539.01577019 7268770.46599761, 2047539.12003661 7268770.46095747, 2047539.22393013 7268770.45080729, 2047539.32720045 7268770.43557153, 2047539.4295988 7268770.41528688, 2047539.53087849 7268770.39000221, 2047539.63079552 7268770.35977844, 2047539.72910919 7268770.32468839, 2047539.82558264 7268770.28481657, 2047539.91998348 7268770.24025906, 2047540.01208426 7268770.19112319, 2047540.10166313 7268770.13752733, 2047540.18850427 7268770.07960061, 2047540.27239848 7268770.01748256, 2047540.35314364 7268769.95132285, 2047540.43054524 7268769.88128085, 2047540.50441681 7268769.80752531, 2047540.57458039 7268769.7302339, 2047540.64086694 7268769.64959283, 2047540.70311677 7268769.56579636, 2047540.76117993 7268769.47904638, 2047540.81491652 7268769.38955187, 2047540.8641971 7268769.29752843, 2047540.90890295 7268769.20319775, 2047540.94892635 7268769.10678709, 2047540.9841709 7268769.0085287, 2047541.01455168 7268768.9086593, 2047541.03999551 7268768.80741948, 2047541.06044109 7268768.70505314, 2047541.07583917 7268768.60180689, 2047541.08615264 7268768.49792945, 2047541.09135667 7268768.39367109, 2047541.09143871 7268768.28928295, 2047541.08639857 7268768.18501653, 2047541.07624839 7268768.08112302, 2047541.06101262 7268767.97785269, 2047541.04072797 7268767.87545434, 2047541.0154433 7268767.77417465, 2047540.98521954 7268767.67425762, 2047540.95012948 7268767.57594395, 2047540.91025767 7268767.4794705, 2047540.86570015 7268767.38506967, 2047540.81656428 7268767.29296888, 2047540.76296842 7268767.20339001, 2047540.7050417 7268767.11654887, 2047540.64292366 7268767.03265466, 2047540.57676394 7268766.9519095, 2047540.50672195 7268766.8745079, 2047540.4329664 7268766.80063633, 2047540.35567499 7268766.73047275, 2047540.27503392 7268766.6641862, 2047540.19123746 7268766.60193637, 2047540.10448748 7268766.54387321, 2047540.01499296 7268766.49013662, 2047539.92296952 7268766.44085604, 2047539.82863885 7268766.3961502, 2047539.73222818 7268766.35612679, 2047539.63396979 7268766.32088224, 2047539.53410039 7268766.29050146, 2047539.43286057 7268766.26505763, 2047539.33049424 7268766.24461205, 2047539.22724798 7268766.22921397, 2047539.12337055 7268766.2189005, 2047539.01911218 7268766.21369648, 2047538.91472405 7268766.21361443, 2047538.81045763 7268766.21865458, 2047538.70656411 7268766.22880476, 2047538.60329378 7268766.24404052, 2047538.50089543 7268766.26432517, 2047538.39961575 7268766.28960984, 2047538.29969871 7268766.3198336, 2047538.20138505 7268766.35492366, 2047538.10491159 7268766.39479548, 2047538.01051076 7268766.43935299, 2047537.91840997 7268766.48848886, 2047537.8288311 7268766.54208472, 2047537.74198996 7268766.60001144, 2047537.65809576 7268766.66212949, 2047537.57735059 7268766.7282892, 2047537.49994899 7268766.79833119, 2047537.42607742 7268766.87208674, 2047537.35591385 7268766.94937815, 2047537.2896273 7268767.03001922, 2047537.22737746 7268767.11381568, 2047537.16931431 7268767.20056567, 2047537.11557771 7268767.29006018, 2047537.06629713 7268767.38208362, 2047537.02159129 7268767.4764143, 2047536.98156788 7268767.57282496, 2047536.94632334 7268767.67108335, 2047536.91594255 7268767.77095275, 2047536.89049872 7268767.87219257, 2047536.87005314 7268767.97455891, 2047536.85465507 7268768.07780516, 2047536.8443416 7268768.18168259, 2047536.83913757 7268768.28594096, 2047536.83905553 7268768.39032909, 2047536.84409567 7268768.49459551, 2047536.85424585 7268768.59848903, 2047536.86948162 7268768.70175936, 2047536.88976627 7268768.80415771, 2047536.91505093 7268768.9054374, 2047536.9452747 7268769.00535443, 2047536.98036476 7268769.10366809, 2047537.02023657 7268769.20014155, 2047537.06479408 7268769.29454238, 2047537.11392995 7268769.38664317, 2047537.16752581 7268769.47622204, 2047537.22545254 7268769.56306318, 2047537.28757058 7268769.64695739, 2047537.35373029 7268769.72770255, 2047537.42377229 7268769.80510415, 2047537.49752783 7268769.87897572, 2047537.57481924 7268769.9491393, 2047537.65546031 7268770.01542585, 2047537.73925678 7268770.07767568, 2047537.82600676 7268770.13573883, 2047537.91550127 7268770.18947543, 2047538.00752471 7268770.23875601, 2047538.10185539 7268770.28346185, 2047538.19826605 7268770.32348526, 2047538.29652444 7268770.35872981, 2047538.39639384 7268770.38911059, 2047538.49763366 7268770.41455442, 2047538.596 7268770.436)

see attached images original circular string is green converted line is red

[claudeunger]

grafic display of the geometries

[claudeunger]

detail of the geometries

[tiiipponen]

Please can somebody hurry this. It makes things very complicated when this kind of essential function is not working.

I'm not sure if problem I describe here is exactly same as original one, but I created following example:

I do in this Python code same things that are in PostGIS source code found from page
​http://postgis.refractions.net/documentation/postgis-doxygen/d6/d41/lwsegmentize_8c-source.html#l00494
in function lwcircle_center:

import math
import sys

p1x=25493284.831099998
p1y=6677883.0782000003
p2x=25493284.841448288
p2y=6677883.0513493409
p3x=25493284.851799998
p3y=6677883.0245000003

temp = p2x*p2x + p2y*p2y
bc = (p1x*p1x + p1y*p1y - temp) / 2.0
cd = (temp - p3x*p3x - p3y*p3y) / 2.0
det = (p1x - p2x)*(p2y - p3y)-(p2x - p3x)*(p1y - p2y)
det = 1.0 / det
cx = (bc*(p2y - p3y)-cd*(p1y - p2y))*det
cy = ((p1x - p2x)*cd-(p2x - p3x)*bc)*det
cr = math.sqrt((cx-p1x)*(cx-p1x)+(cy-p1y)*(cy-p1y))

print("Center({0:f},{1:f}), radius({2:f})".format(cx,cy,cr))

Center(25483200.342228,6673995.406313), radius(10807.909535)

This is clearly incorrect

When I do the same with code found from page
​http://www.mathworks.com/matlabcentral/newsreader/view_thread/294297
made by Roger Stafford

import math
import sys

x1=25493284.831099998
y1=6677883.0782000003
x2=25493284.841448288
y2=6677883.0513493409
x3=25493284.851799998
y3=6677883.0245000003

x21 = x2-x1
y21 = y2-y1
x31 = x3-x1
y31 = y3-y1
h21 = x21**2+y21**2
h31 = x31**2+y31**2
d = 2*(x21*y31-x31*y21)
a = x1+(h21*y31-h31*y21)/d
b = y1-(h21*x31-h31*x21)/d
r = math.sqrt(h21*h31*((x3-x2)**2+(y3-y2)**2))/abs(d)

print("Center({0:f},{1:f}), radius({2:f})".format(a,b,r))

Center(25493495.638481,6677964.308307), radius(225.916095)

This is correct answer.

Can you please hurry repairing this. It is popular function anyway and we want to use it heavily right now.

Thanks

Timo Iipponen
GIS expert
City of Helsinki, Finland

[strk]

Timo, a patch would you the fastest way to get this fixed. You can attach it here. Thank you !

[robe]

By the way, that's the outdated doxygen — please use postgis.net - All listed here: ​http://postgis.net/development/

the 2.1 development one is here: ​http://postgis.net/docs/doxygen/2.1/

[tiiipponen]

Hi

I include here file “lwalgorithm_lw_arc_center.txt”.
It has fixed function “lw_arc_center” from file lwalgorithm.c.
​http://postgis.net/docs/doxygen/2.1/dc/dc0/lwalgorithm_8c_source.html

Original function returned incorrect centroid and radius for arc geometry.
This attached function is minimum correction to fix the problem.

I have tested mathematics of attached function only with Python in Safe Software FME environment.
I also include Python function in file arc_stroking_python.txt for reference.

My Python function has also some extra functionality,
but implementing that in current PostGIS library is more complicated.
I still would like to some day see following functionality in PostGIS ST_CurveToLine function.
Tolerance parameter. In GIS data this is important.
If user gives “Tolerance parameter” value > 0, then “vertex count” parameter means minimum vertices in stroked arc.
For example given tolerance with “vertex count”=0 means, that stroked arc can even be straight line if tolerance allows it.

I am ready to test new fixed ST_CurveToLine function heavily when it is available.

Timo

[tiiipponen]

Includes fixed function “lw_arc_center” from file lwalgorithm.c

[tiiipponen]

Working (FME-) Python function for testing algorithms in file lwalgorithm_lw_arc_center.txt

[tiiipponen]

After thinking a bit more I would like to clarify few comment rows in code:

/* Using cartesian eguations from page ​https://en.wikipedia.org/wiki/Circumscribed_circle */
=>
/* Using cartesian circumcenter eguations from page ​https://en.wikipedia.org/wiki/Circumscribed_circle */

/* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
=>
/* 2 * Cross product, d<0 means clockwise and d>0 counterclockwise sweeping angle */

/* Check colinearity, |Cross product| = 0 */
=>
/* Check colinearity, Cross product = 0 */

[robe]

tiiipponen repackaged as a patch against trunk

[robe]

tiiipponen,

Thanks for the patch, for future reference we require these submissions to be in patch form, rather than just a raw replace of the function. It's easier to just cleanly apply that way to any branch, and also is easier to spot check to see what's changed in it. I repackaged your attachment as a patch as you can see.

I did a quick test against my PostgreSQL 9.5beta2 windows 64-bit instance, and looks like all regress still pass.

pramsey et. al - can you take a look at this and see if it's okay to apply to 2.3, 2.2, 2.1 or we need to add more regression tests for this.

[robe]

P.S. I apologize for my spacing misses. I did my best. I did notice the original code had spaces in some areas (1 at least) where tabs should have been too though so who ever did it in the first place was not perfect

[tiiipponen]

Thanks and sorry. The word patch wasn't familiar to me in this context. Next time…

[robe]

np. You did the hard work .

Thanks, Regina

[strk]

I think such a change should come with an automated testcase clearly showing why the change was needed.

[tiiipponen]

I can try to deliver such a automated testcase if I understand what you want.
I already send two short Python codes.
Would you like to have some PostGIS data dump?
When you look data visually for example with QGIS, you can see bad results.
If you want something else, let me know. I try to give for you everything you need.

[robe]

tiipponen,

Here is our FAQ on how to prepare tests - For lwgeom stuff, we'd need both a cunit test and a postgresql test.

Cunit ​https://trac.osgeo.org/postgis/wiki/DevWikiCUnit

PostgreSQL ​https://trac.osgeo.org/postgis/wiki/DevWikiPGRegress

Ideally we would want the smallest example you can think of that would reproduce the issue (maybe on this ticket, but them side by side so they can compare), so perhaps the example above, but maybe instead of the default number of points to represent an arc, you passed in a smaller number.

Using this form of the function

 ST_CurveToLine(geometry curveGeom, integer segments_per_qtr_circle);

I should add, we normally don't ask this of people reporting bugs or giving patches for fixes, so take it as a BIG compliment, that we think you are capable of doing this .

Normally I'd just do this myself, but I don't work much in curve area and I think pramsey who works on this part is on vacation.

[pramsey]

Committed to trunk at r14425 with a cunit test case.

[pramsey]

Committed to 2.2 at r14426, 2.1 at r14427.

[tiiipponen]

Thank you very much. When I have installed new version of PostGIS and tested that against our data, I let you know the results.

[tiiipponen]

Ok. I have visually checked about 17000 arcs inside different kind of geometries in our data and it looks good to me.

I can also see that this same kind of problem is also in QGIS and in Geoserver. Maybe I have to try to get them fixed too.

But thanks for your help. You are doing a great job.

Timo