PostGIS  2.4.9dev-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 colinear, as determined by equal slopes, then NULL is returned. If the interior point is coincident with either end point, they are taken as colinear.

Definition at line 227 of file lwalgorithm.c.

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

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:328
#define EPSILON_SQLMM
Tolerance used to determine equality.
double y
Definition: liblwgeom.h:328
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
Here is the caller graph for this function: 