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.

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_stroke(), pt_continues_arc(), pta_unstroke(), ptarrayarc_contains_point_partial(), and test_lw_arc_center().

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: