double
lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)	
{
        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;
}