[postgis-tickets] [PostGIS] #3099: ST_curveToLine bad conversion for some arcs

[PostGIS]

#3099: ST_curveToLine bad conversion for some arcs
--------------------------+----------------------------
  Reporter:  claudeunger  |      Owner:  pramsey
      Type:  defect       |     Status:  new
  Priority:  critical     |  Milestone:  PostGIS 2.1.9
 Component:  postgis      |    Version:  2.1.x
Resolution:               |   Keywords:  ST_curveToLine
--------------------------+----------------------------

Comment (by 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[[BR]]
 http://postgis.refractions.net/documentation/postgis-doxygen/d6/d41
 /lwsegmentize_8c-source.html#l00494 [[BR]]
 in function lwcircle_center:

 import math[[BR]]
 import sys[[BR]]

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

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

 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[[BR]]
 http://www.mathworks.com/matlabcentral/newsreader/view_thread/294297
 [[BR]]
 made by Roger Stafford[[BR]]

 import math
 import sys

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

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

 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[[BR]]
 GIS expert[[BR]]
 City of Helsinki, Finland[[BR]]

--
Ticket URL: <https://trac.osgeo.org/postgis/ticket/3099#comment:3>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.