PostGIS  2.1.10dev-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 228 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(), lwcircle_segmentize(), pt_continues_arc(), pta_desegmentize(), ptarrayarc_contains_point_partial(), and test_lw_arc_center().

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

Here is the caller graph for this function: