PostGIS  2.2.7dev-r@@SVN_REVISION@@
 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 213 of file lwalgorithm.c.

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

214 {
215  POINT2D c;
216  double cx, cy, cr;
217  double dx21, dy21, dx31, dy31, h21, h31, d;
218
219  c.x = c.y = 0.0;
220
221  LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
222
223  /* Closed circle */
224  if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
225  fabs(p1->y - p3->y) < EPSILON_SQLMM)
226  {
227  cx = p1->x + (p2->x - p1->x) / 2.0;
228  cy = p1->y + (p2->y - p1->y) / 2.0;
229  c.x = cx;
230  c.y = cy;
231  *result = c;
232  cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
233  return cr;
234  }
235
236  /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
237  dx21 = p2->x - p1->x;
238  dy21 = p2->y - p1->y;
239  dx31 = p3->x - p1->x;
240  dy31 = p3->y - p1->y;
241
242  h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
243  h31 = pow(dx31, 2.0) + pow(dy31, 2.0);
244
245  /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
246  d = 2 * (dx21 * dy31 - dx31 * dy21);
247
248  /* Check colinearity, |Cross product| = 0 */
249  if (fabs(d) < EPSILON_SQLMM)
250  return -1.0;
251
252  /* Calculate centroid coordinates and radius */
253  cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
254  cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
255  c.x = cx;
256  c.y = cy;
257  *result = c;
258  cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));
259
260  LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);
261
262  return cr;
263 }
double x
Definition: liblwgeom.h:312
#define EPSILON_SQLMM
Tolerance used to determine equality.
double y
Definition: liblwgeom.h:312
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55

Here is the caller graph for this function: