PostGIS  2.5.0dev-r@@SVN_REVISION@@

◆ lw_arc_center()

double lw_arc_center ( const POINT2D p1,
const POINT2D p2,
const POINT2D p3,
POINT2D result 
)

Determines the center of the circle defined by the three given points.

In the event the circle is complete, the midpoint of the segment defined by the first and second points is returned. If the points are collinear, as determined by equal slopes, then NULL is returned. If the interior point is coincident with either end point, they are taken as collinear.

Definition at line 227 of file lwalgorithm.c.

References EPSILON_SQLMM, LWDEBUGF, POINT2D::x, and POINT2D::y.

Referenced by lw_arc_calculate_gbox_cartesian_2d(), lw_arc_length(), lw_arc_side(), lw_dist2d_arc_arc(), lw_dist2d_pt_arc(), lw_dist2d_seg_arc(), lwarc_linearize(), pt_continues_arc(), pta_unstroke(), ptarrayarc_contains_point_partial(), and test_lw_arc_center().

228 {
229  POINT2D c;
230  double cx, cy, cr;
231  double dx21, dy21, dx31, dy31, h21, h31, d;
232 
233  c.x = c.y = 0.0;
234 
235  LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
236 
237  /* Closed circle */
238  if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
239  fabs(p1->y - p3->y) < EPSILON_SQLMM)
240  {
241  cx = p1->x + (p2->x - p1->x) / 2.0;
242  cy = p1->y + (p2->y - p1->y) / 2.0;
243  c.x = cx;
244  c.y = cy;
245  *result = c;
246  cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
247  return cr;
248  }
249 
250  /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
251  dx21 = p2->x - p1->x;
252  dy21 = p2->y - p1->y;
253  dx31 = p3->x - p1->x;
254  dy31 = p3->y - p1->y;
255 
256  h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
257  h31 = pow(dx31, 2.0) + pow(dy31, 2.0);
258 
259  /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
260  d = 2 * (dx21 * dy31 - dx31 * dy21);
261 
262  /* Check colinearity, |Cross product| = 0 */
263  if (fabs(d) < EPSILON_SQLMM)
264  return -1.0;
265 
266  /* Calculate centroid coordinates and radius */
267  cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
268  cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
269  c.x = cx;
270  c.y = cy;
271  *result = c;
272  cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));
273 
274  LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);
275 
276  return cr;
277 }
double x
Definition: liblwgeom.h:330
#define EPSILON_SQLMM
Tolerance used to determine equality.
double y
Definition: liblwgeom.h:330
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
Here is the caller graph for this function: