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.

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().

{
        POINT2D c;
        double cx, cy, cr;
        double dx21, dy21, dx31, dy31, h21, h31, d;

        c.x = c.y = 0.0;

        LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);

        /* Closed circle */
        if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
            fabs(p1->y - p3->y) < EPSILON_SQLMM)
        {
                cx = p1->x + (p2->x - p1->x) / 2.0;
                cy = p1->y + (p2->y - p1->y) / 2.0;
                c.x = cx;
                c.y = cy;
                *result = c;
                cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
                return cr;
        }

        /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
        dx21 = p2->x - p1->x;
        dy21 = p2->y - p1->y;
        dx31 = p3->x - p1->x;
        dy31 = p3->y - p1->y;

        h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
        h31 = pow(dx31, 2.0) + pow(dy31, 2.0);

        /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
        d = 2 * (dx21 * dy31 - dx31 * dy21);

        /* Check colinearity, |Cross product| = 0 */
        if (fabs(d) < EPSILON_SQLMM)
                return -1.0;

        /* Calculate centroid coordinates and radius */
        cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
        cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
        c.x = cx;
        c.y = cy;
        *result = c;
        cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));

        LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);

        return cr;
}

Here is the caller graph for this function: